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ABSTRACT 


A  computational  model  has  been  developed  and  applied  for  calculating 
radiant  heat  flux  from  an  electric  arc  using  real  gas  properties  and  taking 
into  account  both  self  absorption  and  space  varying  temperature.  The 
model  requires  knowledge  of  the  temperature  profile  as  an  input.  Tem¬ 
perature  profiles  were  measured  spectrographically  using  an  air  arc  at  one 
and  2. 4  atmospheres.  These  measured  profiles  plus  some  simple  analytic 
ones  which  satisfied  boundary  conditions  were  used  to  calculate  radiant  heat 
flux  at  pressures  from  1  to  300  atmospheres.  The  experimental  program 
also  included  measurements  of  radiant  flux  density  at  2.  2  and  5.  7  atmos¬ 
pheres  for  comparison  with  calculations  based  on  the  measured  profiles. 
Estimates  were  made  of  radiation  losses  in  the  AEDC  20-megawatt  arc  air 
heater  now  under  development. 
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INDEXING  SUBSCRIPTS 

The  following  subscripts  are  used  to  indicate  that  the  subscripted 
quantity  varies  stepwise  with  the  subscript  (listed  in  reverse  order 
of  summation). 
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CL, 

E 

g 

HF 

I 
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Wavenumber 

Index  on  radiating  annulus,  counting  from  the  outermost 

Angle  at  which  absorption  lengths  are  calculated 
(see  Figs.  2.  4  and  2.  5) 

Index  on  absorbing  annulus,  counting  from  the 
outermost  (see  sect.  2.  4) 

GENERAL  SUBSCRIPTS 
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Centerline 

Exterior 

Gas 

High  frequency 
Interior 

Non-black  surface  at  which  net  radiation  is  calculated 
Non-black  surface  radiation  to  surface  i 
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SECTION  I 
INTRODUCTION 


One  of  the  primary  limitations  on  entry  simulation  facilities, 
such  as  hotshot  tunnels,  arc  driven  9hock  tubes  and  electric  arc 
heaters  in  general,  is  radiation  heat  loss.  This  report  describes  a 
project  which  consisted  of  a  computer  study  of  arc  radiation  using 
tabulated  spectral  absorption  properties  of  high  temperature  air, 
combined  with  spectrographic  measurements  of  temperature  profiles 
and  experimental  determination  of  radiation  heat  transfer  rateB. 

The  geometric  configuration  studied  was  an  idealization  of  a 
D„C.  arc  column  in  the  form  of  an  infinitely  long  right  circular 
cylinder  whose  temperature  is  a  known  function  of  radius  only.  The 
basic  approach,  which  makes  use  of  detailed  spectral  absorption  data 
(Refs.  1,  2)  developed  by  Breene  and  Nardone  for  analysis  of  the 
radiating  gas  cap  of  a  re-entry  vehicle,  has  application  to  any 
reasonably  simple  geometry  with  space -varying  temperature. 

An  upper  bound  on  radiation  heat  flux  can  be  calculated  by 
assuming  the  arc  to  be  a  black  body  at  arc  centerline  temperature. 
Another  upper  bound  can  be  established  by  assuming  the  gas  heated 
by  the  arc  to  be  optically  thin  so  that  all  the  energy  radiated  by  each 
infinitesimal  volume  element  reaches  the  chamber  walls.  At  a  given 
pressure,  and  for  a  known  or  assumed  temperature  profile,  radiant 
flux  density  can  then  be  calculated  by  integrating  total  radiance  over 
the  volume  of  a  unit  length  of  arc  column  and  multiplying  by  4  W 
steradians.  The  energy  leaving  this  unit  length  of  arc  column  is,  of 
course,  distributed  over  the  entire  column  length  but,  assuming  an 
infinitely  long  column,  the  energy  radiated  from  a  unit  length  must 
be  equal  to  the  energy  (originating  throughout  the  column)  crossing 
a  unit  length  of  boundary  surface. 

These  upper  bounds  are  perfectly  adequate  as  long  as  radia¬ 
tion  is  only  a  small  part  of  the  total  heat  flux.  However,  at  the 
pressure  and  temperature  levels  which  are  now  the  goals  of,  for 
example,  arc  heater  development,  radiation  is  the  predominant 
heat  transfer  mode  and  a  more  accurate  prediction  is  needed. 

Absorption  and  emission  of  radiant  energy  in  air  depend  on 
temperature,  density  and  frequency.  ThuB,  radiation  energy  from 
a  particular  element  of  volume  will  have  a  spectral  distribution  of 
intensity  which  is  characteristic  of  the  temperature  and  density  of 
that  element.  The  amount  of  energy  reabsorbed  along  a  particular 
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path  will  depend  on  the  path  length  and  also  on  the  absorption  charac¬ 
teristics  (hence  temperature  and  density)  of  each  volume  element 
along  the  path.  A  complete  solution  of  the  radiation  heat  transfer 
problem  with  self -absorbing  gas  therefore  requires  detailed  knowledge 
of  the  radiation  characteristics  of  air  plus  calculation  for  all  volume 
elements  (or  a  representative  sample  if  advantage  can  be  taken  of 
symmetry),  all  possible  paths  for  each  volume  element  and  all  fre¬ 
quencies  for  each  path. 

As  a  practical  matter  some  approximations  are  necessary.  By 
means  of  the  method  described  in  this  report,  sufficiently  tight  upper 
and  lower  bounds  can  be  placed  on  the  effect  of  self -absorption  that  the 
calculation  of  radiant  flux  density  is  limited  in  accuracy  by  the  state  of 
knowledge  of  radiation  properties,  and  of  arc  column  temperature 
profiles. 
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SECTION  II 
ANALYSIS 


2.1  THE  ARC  COLUMN  MODEL 

Consider  an  infinite  cylindrical  column  of  radiating  gas  at 
constant  pressure  with  a  known  or  assumed  radial  temperature  profile. 
Let  the  column  boundaries  be  cold  and  non-reflecting  because,  for 
structural  reasons,  the  wall  must  be  kept  sufficiently  cool  that  re- 
radiation  can  be  neglected  and,  while  cold  compared  with  the  arc 
plasma,  it  is  still  too  warm  to  retain  a  high  reflectivity.  A 
correction  for  the  effect  of  non-black  walls  is  discussed  in 
Section  2,  7. 

Complete  equilibrium  cannot  be  attained  in  the  arc  column, 
which  exists  only  because  of  continuous  energy  transfer  processes. 
However,  a  considerable  body  of  evidence,  both  experimental  and 
theoretical  (Refs.  3,  4,  5,  6  and  7),  indicates  that  local  thermal 
equilibrium  among  electrons,  ions,  atoms  and  molecules  can  be 
assumed,  with  certain  restrictions,  bo  that  equilibrium  properties 
may  be  used  in  the  computation. 

Finkelnburg  (Ref.  3)  has  indicated  that  according  to  existing 
analytical  and  experimental  evidence,  local  thermal  equilibrium 
exists  to  a  good  approximation  in  the  plasma  of  all  arcs  operating  at 
one  atmosphere  or  greater  except  near  the  electrodes  where  strong 
electric  fields  (hundreds  of  volts  per  centimeter)  are  present.  Lochte- 
Holtgreven  (Ref.  4)  concluded  from  his  theoretical  studies  of  the  arc 
column  that,  at  atmospheric  or  higher  pressures,  both  the  temperature 
of  excitation  and  the  "temperature"  of  ionization  are  for  practical 
purposes  equal  to  the  electron  temperature,  with  the  possible  exception 
of  constricted  arcs  having  an  electric  field  intensity  of  300  volts /cm 
or  more.  Cobine  (Ref.  5)  has  shown  that  bulk  gas  temperature  and 
electron  temperature  are  essentially  identical  for  a  constant  current 
discharge  whenever  the  pressure  is  above  17  mm  hg.  Maecker's 
(Ref.  6)  experimental  studies  support  these  conclusions. 

The  arc  heater  used  in  this  study  was  designed  to  isolate 
completely  the  constricted  region  from  electrode  contamination  (a 
major  source  of  column  and  flow  perturbation).  The  only  part  of  the 
arc  column  used  for  this  study  was  the  constricted  region,  well 
separated  from  the  electrodes,  with  a  voltage  gradient  less  than  one 
hundred  volts /cm  (Ref.  8). 
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A  given  temperature  distribution  implies  the  existence  of  a 
steady  state;  each  element  of  volume,  dV,  has  presumably  struck  a 
balance  among  heat  generated,  radiated,  absorbed  and  conducted  in 
and  out.  Each  such  element  within  the  column  Tadiates  in  all  direc¬ 
tions  with  an  intensity  which  is  a  known  function  of  temperature, 
density  and  frequency.  The  assumptions  of  cylindrical  symmetry  and 
constant  pressure  make  temperature  and  density  dependent  on  radius 
only,  so  that  intensity  and  transmissivity  reduce  to  functions  of  radiuB 
and  frequency  only. 

The  arc  column  was  divided  into  a  constant  temperature  core 
and  a  series  of  constant  temperature  annuli*  with  radiant  heat  flux 
originating  in  the  core  and  in  each  annulus  considered  separately  in 
terms  of  equivalent  surface  emissivities.  Radiant  flux  densities  for 
each  surface  are  given  by  an  equation  of  the  form 

q/A  =  |  c  WQ  dn  (2.  1) 

J  n  B,n 


where  W3  n  is  the  monochromatic  emissive  power  of  a  black  body  at 
the  radiator  temperature  and  n  is  the  wavenumber.  Over  a  small 
but  finite  wavenumber  interval  the  flux  density  is  of  the  form 

(q/A)  =  a  T4  e  T)  (2.2) 

tin  n  n 

where  Cn  is  the  average  emissivity  within  the  wavenumber  interval  and 
7Jn  is  the  fraction  of  the  total  energy  radiated  by  a  black  body  at  tem¬ 
perature  T  which  is  radiated  within  the  interval  An. 


Each  annulus  was  denoted  by  the  index  number  j  assigned  to  its 
exterior  surface,  counting  from  the  outside  surface  of  the  arc 
column  inward  to  the  core. 

Wavenumber  is  the  reciprocal  of  wavelength  and  is  therefore  pro¬ 
portional  to  frequency.  Tabulation  of  properties  by  wavenumber 
proves  convenient  so  the  term  will  be  used  interchangeably  with 
frequency  in  the  sense  of  frequency  dependence. 
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Upper  and  lower  limits  were  placed  on  the  fraction  of 
emitted  radiant  energy  which  actually  leaves  the  arc  column  by 
considering  two  assumed  sets  of  absorbing  paths.  One  set  of  paths 
was  so  chosen  that  absorption  is  always  less  than  or  equal  to  the  true 
absorption  while  along  the  other  set  of  paths  radiation  absorption 
was  always  greater  than  or  equal  to  the  true  value. 

2.2  INFINITE  HALF  CYLINDER  ABSORPTIVITY 

The  basic  element  from  which  the  absorption  model  was 
built  up  is  an  infinite  half -cylinder  of  gas,  of  radius  R,  absorbing 
radiation  from  an  element  of  area  on  its  axis.  Associated  with  the 
cylinder  is  a  linear  absorption  coefficient  k  which  is  a  function  of 
radial  distance  only.  The  angles  shown  in  Fig.  2.1  (a)  are  related  by 


cos  v  =  sin  Q  cos  1 1> 


(2.3) 


From  Fig.  2.1(b),  the  solid  angle  dw  can  be  written 


dto  a  d$d  0 

Xj 


(2.4) 


If  the  surface  element  dA  is  radiating  with  an  intensity  I  (energy 
flux  per  unit  solid  angle)  and  absorption  takes  place  along  the  path 
according  to  Beer's  law,  then 


d  q  =  (I  cos  ydA)  dw  exp 


(2.5) 


Integration  of  Eq.  (2.  5)  over  the  half  cylinder  would  yield  the  total 
radiant  flux  originating  from  dA  and  escaping  the  half  cylinder. 
With  the  further  definition 


(kR) 


*  R 

kdr 

o 


I 

sin  Q 


f  L 

kd£ 


(2.6) 


the  radiant  flux  density  from  dA  which  escapes  the  half  cylinder  is 
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Figure  2.  i  Absorption  by  an  Infinitely  Long  Half  Cylinder. 
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dq 

dA 


=  41 


r  ff/2 


cob  0  d  0 


viz 


KR 
sin  0 


sin  & 

(2.7) 


For  the  case  of  no  absorption 
q/A  =  ir  I 


(2.8) 


Therefore  if  kR  j^f(  0)  transmissivity  and  absorptivity  are  given  by 


4 

t  =  1  -  a  =  — 

*  IT  fZ 

COS  0d0 

•  V  /  2 

exp 

r  (kR)n  1 

n  n  V 

• 

o  • 

o 

sin  9  1 

,  Znjn 
sin  9d0 


(2.9) 


■where  subscript  n  denotes  dependence  on  wavenumber.  Equation 
(2.  9)  cannot  be  integrated  inclosed  form  but  numerical  evaluation 
is  straightforward  and  the  r£sult  is  Bhown  in  Fig.  2.2,  Note  that 
an  and  T  depend  only  on  (kR^  hereafter  called  the  absorption 
length.  $hen  kn  is  assumed  to  undergo  a  serieB  of  Btep  changes 
with  radius,  Eq.  (2.  6)  is  replaced  by  an  appropriate  summation 
but  Eq,  (2,9)  is  unchanged. 

2.3  ABSORPTION  OF  EXTERIOR  RADIATION 

An  element  of  the  exterior  surface  of  any  annulus  but  the 
outermost  "sees"  the  annuli  exterior  to  it  as  a  segment  of  a  cylinder. 
The  segment  was  approximated  by  two  or  more  sectors  of  cylinders 
centered  at  the  radiating  surface  element  is  shown  in  Fig.  2,3,  In 
this  way,  calculation  of  absorption  became  simply  a  matter  of  de¬ 
termination  of  absorption  lengths  plus  application  of  Eq.  (2.9). 

The  boundaries  of  a  typical  sector  are  at  angles  0and  0+/ 1  0 
from  the  normal  to  the  radiating  surface  element.  It  is  clear  from 
Fig.  2,3  that  an  upper  limit  transmission  approximation  --  all 
absorption  lengths  less  than  or  equal  to  the  true  absorption  lengths 
--  corresponds  to  the  case  where  sector  coincides  with  segment  at 
the  0boundary.  Likewise,  a  lower  limit  on  transmissivity  corre¬ 
sponds  to  sector -segment  coincidence  at  the  0  +  £0  boundary. 
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Figure  2.  2  Transmissivity  of  an  Infinitely  Long  Half  Cylinder. 
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LOWER  LIMIT  ON  UPPER  LIMIT  ON 


Figure  2.  3  Absorption  of  Exterior  Radiation.  The  j  =  3  Annulus  is 
Shown  as  a  Typical  Example. 
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The  absorption  geometry  calculation  must  be  done  separately 
for  all  but  the  outermost  (j  =  1)  annulus.  Referring  to  Fig.  2,4 
there  are  q  =  j  -  1  sector  radii*  R.  ^  exterior  to  the  jth  annulus 
within  each  sector.  Calculation  of  each  Rj^  q  is  straightforward 
trigonometry  since  r;  and  r  are  known  and  4  is  equal  to  either  ip  for 
the  upper  limit  calculation  Ip  +  A^  for  the  lower  limit  calculations. 
Absorption  lengths  are  then 


j*  4 


'’J.  q.  +  1 


(2.  10) 


The  magnitude  of  A$  depends  on  the  number  of  sectors  into 
which  the  quadrant  0  £  Ip  £  was  divided.  Let  s  be  the  number  of 
lectors,  so  that 


A  ip  = 


TT 

2s 


(2.11) 


Transmissivity  is  then  given  by 


fl/2 


n.J 


®  [  |  r  r-(kR)n  •  <k 

= -^T  |sin(i  A  0>-sin  [(i-i)  A  exp  si 

i=l 


(2.12) 


Subscript  j  on  sector  radius,  absorption  length  and  transmissivity 
refers  to  the  respective  quantities  as  "seen"  from  the  annulus.^ 
Subscript  q  is  used  to  denote  annuli  absorbing  radiation  from  the  j 
annulus,  again  counting  inward  from  the  outermost  radius.  Capital 
R  is  used  to  distinquish  sector  radii  from  annulus  radii  and  the  angle 
A  ip  designates  the  angle  at  which  absorption  lengths  for 
a  particular  sector  were  calculated. 
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Upper  and  lower  limits  on  transmissivity  depend  on  choice  of  0 


i  =  i  A  0;  T=  (T  .) 

j  LL 

=  (i-i)A  T=  (r  ) 

n*J  E, UL 


{2.13a) 


where  subscripts  E,  LL  and  UL  refer  to  exterior,  lower  limit  and 
upper  limit  respectively. 

2.4  ABSORPTION  OF  INTERIOR  RADIATION 


The  parts  of  annuli  "visible"  to  a  radiating  element  of  in¬ 
terior  surface  are  no  longer  segments,  but  similar  sector  approx¬ 
imations  can  be  applied.  Calculation  of  absorption  lengths  is  some¬ 
what  complicated  by  the  fact  that,  depending  on  0 ,  a  radiation  path 
can  cross  an  annulus  boundary  twice,  once,  or  not  at  all,  corre¬ 
sponding  to  the  existence  of  two,  one,  or  no  solutions  to  the  triangle 
specified  by  r  .  +  ^  r  and  0,  Fig.  2.  5.  Sector  radii*  are  given  by 


Rj.  0.q 


sin  0 


’ 

- 1  /  r  \  • 

sin 

0+  sin  (  -■  -  sin  0) 

V  rq  J\ 

(2.14  a) 


when 


r  sxn  0  s  r 

j*1  q 


and  by 


R 


j,  0,  q  sin  0 


|  sin  ^  0  +  v  -  sin  ^  ^  —  sin  ^ 


{2.  14b) 


when 


r  sin  0  < 
J+l 


s  r 


j  +  1 


Sector  radii  equal  to  or  larger  than  r  .  are  designated  R  while  those 
smaller  than  r  are  designated  R".  J 
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Figure  2.  5  Geometry  of  Absorption  of  Interior  Radiation  Illustrated 

for  Radiation  from  the  Outermost  Annulus  (j=l)  with  q=j+2. 
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Defining  q1  as  the  maximum  value  of  the  index  q  {corre¬ 
sponding  to  the  smallest  r  )  for  which  two  solutions  exist'  the 
absorption  lengths  are  givln  by 


(kR) 


n,  j, 


-1 


(2.15) 


+  k 


n,  q’ 


) 


q' 

*  2 

q  =  j+2 


k 

n,  q- 


1 


( 


R.  JL 
j.  0.  q 


Because  absorption  coefficient  is  not  necesarily  a 
monotonically  increasing  function  of  temperature,  the  extremals  of 
of  absorption  length  within  any  sector  do  not  always  correspond  to 
the  values  at  sector  boundaries  Ip  and  0  +  Alp.  Therefore,  a  third 
computation  of  was  made  within  each  sector  at 

4>  =  ip  +  A  ip/2,  The  smallest  and  largest  of  these  were  used  for 
computing  interior  transmissivities  (Tn>  j)  l( 
respectively  from  Eq.  (2.12).  ’  'LL 


“d(rn,j)i. 


UL 


2.5  EMISSIVITY  CALCULATIONS 

Radiation  from  the  exterior  and  interior  surfaces  of  each 
constant  temperature  annulus  was  also  considered  separately. 
Emissivity  of  the  exterior  surface  was  obtained  by  modifying  the 
limits  on  the  absorptivity  integral  for  a  complete  cylinder.  At  the 
interior  surface  emissivity  was  obtained  using  a  modification  of 
the  transmissivity  calculation  already  discussed. 
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The  basis  for  calculating  emissivity  of  the  exterior  surface 
is  the  integral  expressing  absorptivity  of  an  infinitely  long  gaseous 
cylinder  relative  to  an  element  of  its  surface.  The  integral  is 
given  by  Jakob  (Ref.  9)as 


CL 


+  77  /2  +  ?r/2 

p  /* 


/■* 

/* 

_  L 

2 

f  -kr 

COS  Xcos  0 

l  e  dr 

•/  J_  J 

-  77/2 

-  tt/2 

_  J  o 

dXd0 

(2.16) 


where  0  is  defined  in  Fig.  2.  6  and  Xis  the  angle  between  the  vertical 
and  L  measured  in  a  plane  passing  through  the  surface  element 
normal  to  the  plane  of  Fig.  2.  6. 

Instead  of  evaluating  the  complete  double  integral,  a  closed 
form  solution  was  obtained  for  the  case  where  kL,  is  sufficiently 
small  that  the  inner  integral  reduces  to 


kL  <<  1 


(2.17) 


This  solution  was  extended  to  optically  denser  conditions  using  the 
mean  beam  length  concept  (Ref.  10). 


Replacing  the  cylinder  by  a  cylindrical  annulus.  Fig.  2,6, 
only  the  unshaded  part  of  the  annulus  is  considered  because  radiation 
which  passes  through  the  hollow  is  included  in  the  emissivity  of 
the  interior  surface.  The  integration  of  Eq.  (2.  16)  over  0must 
therefore  be  performed  in  two  steps: 


0  =£050  :  L  = 


± 


cos  X 


cos  0  -  \  .  2_  .  2 

y  sin  p  -  sin  0 


(2.  18a) 


COS  0 


(2. 18b) 

where  S  is  the  angle  at  which  L  is  tangent  to  the  inner  surface  and  is 
related  to  the  ratio  of  annulus  radii  by 


z  ib  z  — 
v  2 


2r . 


L  = 
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Figure  2.  6  Absorption  Paths  for  Calculation  of  Exterior  Emissivity  from 

an  Infinitely  Long  Cylindrical  Annulus.  Only  the  Unshaded  Part 
of  the  Annulus  is  Directly  Visible  to  the  Surface  Element. 
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sin  jS 


(2.19) 


Combining  Eqs.  2.  16,  2. 17  and  2. 18  the  absorptivity  of  an  annulus 
can  be  expressed  by 


2kr . 


ttj  = 


+  JL 
r  z 


IT 

2 


cos  Xd  X 


(2.20) 


- 

"  ^  /  j  2  2  \ 

r  IT/  2 

(  cos  if)  -  \ /sin  /S'-  sin  </)  J  cos  tb  d  lb  +  2 

2 

cos  lb  d  0 

J 

o  \  '  / 

.8 

Integrating  and  dividing  by  2krt,  which  is  the  limiting  absorptivity 
of  a  cylinder  of  radius  r.  J 


a 


J 


2kr , 
J 


•  2  A 
sin  p 

2 


0  sin  2  ft 
IT  '  2ir 


(2.21) 


Equation  (2.21)  is  plotted  in  Fig.  2.  7.  If  the  factor  2  is 
changed  to  1.9*  the  mean  beam  length  of  a  cylinder  is  obtained  and 
the  absorptivity  over  the  range  fr.om  optically  thin  to  black  body  is 
given  within  +5%  (Ref.  10)  by 


a 


[  1  -  exp  (-1.9  kr.)| 
J 


(2.22) 


Combining  Eqs.  (2.21)  and  (2.  22)  and  noting  that  the  assumption  of 

local  equilibrium  implies  a  =  e ,  the  following  expression  for  "exterior" 

emissivity  (  €  .  1  results 

J»  n  E 
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Figure  2.  7 


Absorptivity  of  a  Cylindrical  Annulu 
of  a  Cylinder  (both  infinitely  long). 


s  Relative  to  Absorptivity 
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1  -  exp 


1.9k  r 

n  J 


2a 
Bin  p 

2 


sin23  \ 

2  it  j 


(2.23) 


Equation  (2.23)  includes  the  emissivity  of  the  cylindrical  core  as  the 
special  case  sin  #=  0. 


The  interior  emissivity  of  an  annulus  is  identical  with  the  ab¬ 
sorptivity  of  a  single  annulus  relative  to  a  radiating  element  at  its 
inner  surface.  The  sum  of  absorptivity  and  transmissivity  is  unity 
so  transmissivity  of  a  single  annulus  can  be  calculated  from  Eq.  (2.  12) 
using  absorption  lengths  given  by 


<kR>=(kn,j>  IV.rf.j’ 


(2.24) 


The  arithmetic  mean  of  the  upper  and  lower  limits  on  rwas  used  to 
calculate  interior  emissivity 


1 


T  -t-  T 
1  UL  'LL 


2 


(2.25) 


2.6  COMPLETE  EXPRESSION  FOR  RADIANT  HEAT  FLUX 


Radiant  heat  flux,  expressed  as  power  per  unit  surface  area  of 
arc  column  boundary,  may  now  be  written  as  a  double  summation  of 
terms  like  Eq.  (2.2),  each  multiplied  by  the  appropriate  transmissivity 


q/A  =  a 


I 

n 


m 


j=l 


r. 


(  T  .) 

(«  : ) 

r  m-1 


y  T 4p4i-)U  ,)  (r  .)  («  .) 

L  2L  J  x  O  '  x  rn,j'  n,j'  x  n,j' 

j=l 


(2.26) 


19 


AEDC-TR-65-U 


The  radius  of  the  arc  column  is  6,  m  is  the  number  of  constant 
temperature  subdivisions  {core  plus  annuli)  and  T  and  €  are  given  by 
Eqs.  (2.  13),  (2.23)  and  {2.25),  Black  body  energy  fraction,  t},  i9 
propo rational  to  the  integral  of  the  Planck  radiation  function  over  each 
wavelength  interval  An.  Methods  used  for  calculation  of  t]  are  dis¬ 
cussed  in  Appendix  A. 

2.7  CONSTRICTED  ARC  WITH  NON  BLACK  WALLS 

As  pointed  out  in  Section  2.  1  the  arc  chamber  wall  will  have  a 
low  reflectivity  {high  emissivity).  However,  wall  emiseivity  will  not 
be  unity*  as  has  been  assumed  thus  far,  and  a  means  of  evaluating  its 
effect  on  the  magnitude  of  the  radiant  heat  flux  actually  absorbed  by 
the  wall  is  desirable. 

A  complete  calculation  of  the  radiant  interchange  between  the 
arc  and  the  constrictor  walls  adds  a  very  substantial  complication  to 
the  analysis  already  presented  and,  in  addition  requires  data  on  the 
frequency  and  angular  dependence  of  wall  surface  emissivity. 
Fortunately  the  effect  of  a  non -black  wall  can  be  accounted  for  quite 
closely  in  the  following  way. 

Radiant  interchange  between  a  non -black  wail  and  a  radiating 
gas  can  be  expressed,  following  Eckert  (Ref.  10) 

(q/A).  =  -  B.  +  Z  Fi.k<'  -  “i  -k,g>  Bk  +  <*/A>g  (2-27> 

k 

where  (q/A).  is  the  net  radiant  flux  density  absorbed  at  surface  i, 
and  B^  are  the  total  radiant  flux  densities  leaving  surface  i  and 
surface  k  respectively,  is  a  view  factor  between  surface  i  and 

each  surface  k,  ai_kf  g  is  the  absorptivity  of  the  gas  between  sur¬ 
face  i  and  each  surface  k  as  determined  by  integration  over  the  solid 
angle  corresponding  to  the  particular  F.^,  and  (q/A)  is  the  radiant 
flux  density  from  the  gas  incident  on  surface  i.  Equation  2.27 
expresses  the  net  radiant  flux  density  to  surface  i  as  the  algebraic 
sum  of  radiant  flux  density  leaving  surface  i,  radiant  flux  density 
arriving  from  surfaces  k  and  radiant  flux  density  from  the  gas. 


McAdams  (Ref.  17)  gives  a  total  normal  emissivity  of  0.  57  for  a 
copper  "plate  heated  at  1110^". 
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Let  surface  i  be  an  infinitesimal  element  of  the  chamber  wall 
and  consider  as  a  single  surface  k  the  remainder  of  the  infinitely 
long  cylindrical  arc  chamber  wall.  Then 


(q/A)  =  (q/A)  ~  (2.28) 

8  c 

where  (q/A)a  is  the  radiant  flux  density  from  the  arc  as  given  by 
Eq.  2.  26  and  R  is  the  constrictor  radius.  Also 

C 

B.  =  B,  =  B  (2.29) 

i  k  w 

F.  -  k  =  1  (2.  30) 


Z  Fi-k<1 

k 


a  ) 

i-k,  g 


B,  =  (1 


a  )  B 
w.  g  w 


(2.31) 


where  B  is  the  radiant  flux  density  leaving  the  wall  at  any  point  and 

Qfw  is  the  "average"  absorptivity  of  the  gas  between  surface  element 

i  and  all  other  parts  of  the  wall  surface.  Further,  defining  B7  as  the 

total  radiant  flux  density  incident  on  the  wall  at  any  point  we  can  write, 

where  f  is  the  emissivity  of  the  wall 
w 


B  =  (  CT4  +  (1  -  f  )  B' 

WWW  w  w 


(2.32) 


By  definition 


(q/A) 


w 


B'  -  B 
w  w 


(2.33) 


so  that,  neglecting  radiation  at  wall  temperature. 


(q/A) 


w 


Combining  eqs.  2,27,  2.28,  2,31  and  2,  34  we  obtain 


(q/A) 

w 

(q/A)  ( 

a 


(2.34) 


(2.35) 
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which  has  been  plotted  as  Fig.  2.  8. 

The  quantity  is  an  '‘average",  over  all  frequencies  and 

angles,  of  the  absorptivity  of  the  gas  relative  to  radiation  which  has 
been  reflected  one  or  more  times  from  the  constrictor  wall.  Reflec¬ 
tion  may  itself  vary  with  frequency  and  angle,  and  multiple  relfection 
will  be  more  likely  at  those  frequencies  at  which  gas  absorption  is 
relatively  low.  However,  if  the  arc  fills  the  constrictor,  aw,g  cannot 
differ  greatly  from  the  average  emiBsivity  of  the  arc  column  and,  as 
arc  radius  becomes  small  compared  to  wall  radius,  awFg  approaches 
zero.  A  reasonable  estimate  of  aw,g  and  Fig.  2.8  shows  that  the 
correction  is  relatively  insensitive  to  the  precise  value  of  aw,g* 
particularly  for  the  relatively  high  values  of  emissivity  expected  in 
practice. 

All  results  in  Sections  3,  4,  and  5  are  presented  in  terms  of  heat 
flux  from  the  arc  and  a  wall  emissivity  correction  may  be  applied  as 
required. 
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Figure  2.  8  Influence  of  Arc  Chamber  Emissivity  on  Radiant  Heat 
Transfer, 
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SECTION  111 

DIGITAL  COMPUTATION 


3.1  PROGRAM  ARCRAD 

The  FORTRAN  computer  program  which  enabled  an  IBM  7090 
digital  computer  to  produce  solutions  to  Eq.  2.23  is  called  ARCRAD 
and  is  included  as  Appendix  B.  It  requires  as  input  information  the 
temperature,  density  and  inner  and  outer  radii  of  each  annulus,  ab¬ 
sorption  coefficient  data  and  some  programming  data  discussed  in 
Appendix  B.  Its  standard  output  includes 

1,  A  repeat  of  the  temperature-density-radius  input, 

2,  Radiant  flux  density  based  on  upper  and  lower  limit 
transmissivities  tabulated  by  wavenumber  interval  and 
also  summed  over  all  wavenumbers, 

3.  Radiant  flux  density  with  all  transmissivities  set’ equal 
to  unity,  again  both  tabulated  by  wavenumber  interval 
and  summed. 

4.  Optionally,  an  output  tape  which  can  be  converted 
automatically  to  a  graph  of  tabulations  2)  and  3)  above. 

3.2  LIMITS  ON  FREQUENCY  RANGE  OF  CALCULATIONS 

Although  radiant  energy  is  distributed  throughout  the  elec¬ 
tromagnetic  spectrum,  all  but  a  negligible  fraction  is  contained  with¬ 
in  rather  narrow  limits  which  are  a  function  of  temperature.  These 
limits  are  conveniently  expressed  in  terms  of  a  parameter  v  defined  as 


v 


°2n 


(3.1) 


where  C2  is  Planck’s  second  radiation  constant  and  is  equal  to 
1.4380  cm-  K,  For  a  black  body  radiating  at  a  given  temperature, 
only  0,25%  of  the  total  energy  is  radiated  at  frequencies  corresponding 
to  v  <  0.  3  or  to  v  >  13,.  Therefore  the  only  energy  considered  was 
that  at  frequencies  (wavenumbers)  bracketed  by  the  above  values  of  v. 
Figure  3.  1  shows  these  limits,  along  with  the  curve  v  =  1,  which  was 
the  dividing  line  for  application  of  two  different  methods  for  calcula¬ 
ting  rj  (see  Appendix  A). 
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The  solid  outline  of  Fig.  3.  1  represents  the  range  of  temper¬ 
ature  and  frequency  of  the  Breene  and  Nardone  results  (Refs.  1,  2,1 
plus  additional  as  yet  unreported  work.  Data  above  61,  500cm-*are 
not  precise  but  are  known  to  provide  an  upper  bound  within  a  factor  of 
2.  Temperatures  near  the  centerline  on  the  measured  profiles  ex¬ 
ceeded  the  15000°K  upper  limit  on  available  high  frequency  data.  An 
upper  bound  on  the  possible  contribution  in  the  high  frequency  range 
was  therefore  established  by  considering  annuli  at  temperatures  above 
15000°K  to  be  black  bodies  in  the  wavenumber  range  for  which  data 
was  not  available. 

m  /  r  \ 

(q/A)HF  °  j  T  j  (  6^  )  (D6l,  5000*j  (3.2) 

3=1 

J  J16, 000 

where  D  is  defined  by  Eq.  A.  2  and  is  here  the  fraction  of  black  body 
energy  above  61,  5000  cm-*.  The  dependence  of  D^  ^  on  temper¬ 
ature  is  shown  in  Fig.  3.2,  Both  D  and  (r/6)  are  generally  small  in 
this  case  so  use  of  the  upper  bound  does  not  greatly  increase  the 
'uncertainty  of  the  results. 

3.3  PRELIMINARY  PROCESSING  OF  ABSORPTION 
COEFFICIENT  DATA 

The  Breene  and  Nardone  data  from  1300  cm  1  to  62,000  cm  * 
is  recorded  on  magnetic  tape  but  the  tape  contains  much  information 
extraneous  to  this  study  and  the  absorption  coefficients  are  tabulated 
at  unnecessarily  fine  frequency  intervals.  A  preliminary  processing 
of  these  tapes  was  therefore  necessary  to  obtain  a  data  tape  containing 
linear  absorption  coefficients  tabulated  at  500  cm-*  intervals. 

Data  above  62,  000  cm  *was  available  only  in  printed  form 
so  punched  cards  were  prepared  which  were  then  converted  to  mag¬ 
netic  tape.  The  wavenumber  interval  was  increased  to  1000  cm-* 
because  at  higher  frequencies  a  given  increment  represents  a  smaller 
fraction  of  the  total  energy. 

3.4  TEMPERATURE  AND  DENSITY  PROFILES 
USED  IN  COMPUTATION 

As  has  been  stated,  temperature  distribution  in  the  arc  column 
is  presumed  known.  With  the  assumption  of  local  thermal  equilibrium 
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TEMPERATURE  °K 


Figure  3.  2  Fraction  of  BLack  Body  Energy  Emitted  at  Wave  Numbers 
Greater  Than  61, 500  cm"1. 
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and  a  given  pressure  level,  the  density  distribution  can  then  be  deter¬ 
mined  from  a  table  or  chart  of  thermodynamic  properties,  A  limited 
amount  of  temperature  profile  data  was  available  from  direct  spec- 
troscoptic  measurement  (see  Section  5,2  for  further  discussion)  which 
provided  a  basis  for  comparison  of  results  with  direct  measurements 
of  radiant  heat  flux. 


Based  on  the  results  of  direct  measurements  plus  application 
of  reasonable  boundary  conditions,  a  series  of  simple  analytic  (poly¬ 
nomial)  profiles  were  also  used  for  computation.  These  profiles 
were  chosen  to  provide  a  systematic  comparison  of  the  effect  on 
radiant  heat  flux  of  temperature  level  and  fullness  of  the  temperature 
profile.  These  analytic  profiles,  shown  in  Fig,  3,3  along  with  a  pro¬ 
file*  measured  by  Maecker,  (Ref.  6)  were  of  the  form 


ref 


T  -  T  , 
CL  ref 


(3.3) 


The  spectral  absorption  coefficient  data  were  tabulated  at 
temperature  intervals  of  1000°K  and  decade  density  intervals.  In 
order  to  avoid  double  interpolation  for  temperature  and  density,  the 
thickness  of  each  annulus  was  selected  in  such  a  way  that  the  weighted 
avarage  temperature  of  each  annulus  was  a  tabulated  temperature, 
see  Appendix  C.  Absorption  coefficients  were  then  interpolated 
logarthmically  with  density.  Calculation  of  weighted  average  annulus 
temperatures  requires  a  polynomial  representation  of  the  temperature 
profiles  so  measured  profiles  were  reduced  to  this  form  using  a 
modification  of  a  standard  least  squares  polynomial  curve  fit  program, 
Appendix  D, 

3.5  TRANSMISSIVITY  APPROXIMATION 

Numerical  evaluation  of  the  integral  appearing  in  Eq.  2,  12 
leads  to  the  curve  of  Fig.  2.2  with  (kR)n  j  ^  as  a  parameter.  Rather 
than  perform  the  numerical  integration  for  each  different  (kR)n  j  ^ 
empirical  approximations  were  used  in  Program, ARCRAD  which  took 
advantage  of  the  similarity  of  Fig.  2.2  to  a  simple  exponential  decay. 
The  approximations  were  made  over  three  regions  as  follows; 

1)  0  *  kR  *  1.1  T  =  exp  -1.206  (kR) 
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2)  l.l<kRS16  T  =  t  exp  C-1.206  {kR)  ]}  Cl.  0774-.  09539  (kRl 

+  .055577  (kR)2  -  .004291  (ER j3 

+  .00020826  (Kl)4} 

3)  kR  >  16  T  =  0 

Introduction  of  this  approximation  resulted  in  roughly  halving  the  com¬ 
puter  time  per  case. 
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SECTION  IV 
COMPUTER  RESULTS 

4.1  WAVENUMBER  RANGE  OF  RESULTS 

Because  of  the  limitations  on  available  data  discussed  in 
Section  3.2,  curves  presenting  results  of  calculation  are  labeled  as 
follows  to  indicate  the  wavenumber  range  they  cover 


Label 

Meaning 

W  a  venumber 
Range 
(cm-!) 

T  emperature 
Range 
(°  K) 

LF 

Low  Frequency 

up  to  61, 5000 

3,  000  -  25,000 

HF 

High  Frequency 

61, 500  and  up 

9,  000  -  15,000 

BBL 

Black  Body  Limit 
given  by  Eq,  3,  2 

61, 500  and  up 

16,000  -  25,  000 

Table  4.  1  -  Frequency  Range  of  Results 

As  will  be  discussed  in  Section  5,4  the  break  at  61,  500  cm 
is  also  convenient  from  the  standpoint  of  comparisons  with  experi¬ 
mental  measurements  since  this  frequency  is  just  about  the  cutoff 
frequency  of  the  fused  quartz  window  used  to  protect  the  radiant  heat 
flux  gage. 

4.2  CONVERGENCE  OF  UPPER  AND  LOWER  BOUNDS 

The  computational  model  described  in  Section  2  of  this  re¬ 
port  was  designed  to  bracket  the  effect  of  self  absorption  on  radiant 
heat  flux  from  an  arc  column.  The  effectiveness  of  this  approach 
depended  on  just  how  close  to  each  other  were  the  upper  and  lower 
bounds  and  this  could  not  be  established  without  actually  performing 
some  numerical  calculations. 

As  originally  conceived  and  proposed,  the  absorption  model 
postulated  a  lower  limit  on  absorption  (maximum  radiant  heat  flux) 
in  which  all  radiation  was  considered  to  be  "exterior"  in  the  sense  of 
section  2.3.  For  this  case  the  absorption  paths  were  represented  by 
half -cylinder s  whose  radii  were  equal  to  the  difference  between  annulus 
outside  radii  and  the  arc  column  boundary  radius.  The  upper  limit  on 
absorption  (minimum  radiant  heat  flux)  was  established  by  assuming 
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all  radiation  to  be  "interior"  and  absorbed  by  half -cylinders  whose 
radii  were  equal  to  the  sums  of  the  annulus  radii  and  the  arc  column 
boundary  radius*  Results  with  this  model  gave  an  upper  bound  prac¬ 
tically  indistinguishable  from  the  optically  thin  case  and  a  lower  bound 
which  was  lower  by  an  order  of  magnitude. 

A  second  computational  model  was  then  developed,  as  des¬ 
cribed  in  detail  in  Section  2.  Division  of  the  half  cylinders  into  sec¬ 
tors  and  separate  calculation  of  "exterior"  and  "interior"  radiation 
for  both  upper  and  lower  bounds  greatly  reduced  the  spread  between 
them.  By  retaining  the  number  of  sectors  as  variable,  more  sectors 
could  be  used  at  higher  pressures,  where  the  spread  iB  large,  than  are 
needed  at  lower  pressures. 

The  additional  flexibility  achieved  by  retaining  the  number 
of  sectors  as  a  variable  can  be  used  in  two  ways.  Figure  4. 1  illus¬ 
trates  the  effect  of  varying  the  number  of  sectors  while  keeping  all 
else  constant.  Within  the  limits  of  roundoff  error,  increasing  the 
number  of  sectors  improves  the  approximation.  A  price  is  paid  of 
course,  in  additional  computer  time.  The  form  of  convergence  in¬ 
dicates  that  an  average  of  upper  and  lower  bound  calculations  will 
give  an  excellent  approximation  to  the  "exact"  result. 

In  addition  to  justifying  an  average  of  the  two  limit  calculations 
as  representative  of  the  "exact"  result.  Fig.  4.  1  is  a  useful  guide  to 
economic  choice  of  the  number  of  sectors.  At  pressures  equal  to  or 
greater  than  100  atm,s  was  set  equal  to- 3.  Below  100  atm  a  value  of 
s=2  sufficed. 

4.3  RESULTS  WITH  ASSUMED  TEMPERATURE  PROFILES 

If  one  starts  with  the  postulate  that  the  arc  column  is  axially 
symmetric  and  that  the  temperature  decreases  monotonically  from 
a  peak  at  the  centerline  to  some  relatively  low  temperature  at  the 
outermost  radius,  than  the  general  shape  of  the  profile  is  reasonably 
well  determined.  This  is  true  even  though  local  bumps  may  exist 
due  to  changes  in  ionization  level  or  other  phenomena. 

Let  us  define  a  full  profile  as  one  in  which  arc  temperature 
remains  near  arc  centerline  temperature  until  very  near  the  outer 
boundary.  Fullness  is  dependent  on  details  of  the  local  energy  bal¬ 
ance  and  is  therefore  certain  to  change  with,  for  example,  changes  in 
the  relative  importance  of  the  conduction  and  radiation  modes  of  heat 
loss  to  the  boundary. 
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SECTORS  PER  QUADRANT  -  S 


Figure  4.  1  Convergence  of  Upper  and  Lower  Limit  Calculations  as 
Number  of  Sectors  Increases. 
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The  family  of  curves  given  by  Eq.  3.3  satisfies  the  afore¬ 
mentioned  postulates  and  contains  a  single  parameter,  the  exponent  x, 
to  characterize  the  fullness.  The  effect  of  fullness  may  then  be  judged 
by  a  plot  such  as  Fig.  4.2  which  shows  radiant  flux  density  vs.  x  at  3 
pressure  levels.  Fig,  4.2  also  shows,  for  comparison,  the  asymptotic 
value  which  each  curve  approaches  as  x  becomes  very  large.  The 
asymptotes  were  computed  assuming  a  cylinder  of  gas  uniformly  at 
centerline  temperature.  The  limiting  case  of  a  blck  body  at  center - 
line  temperature  is  also  included.  Lacking  knowledge  of  temperature 
profiles  at  the  higher  pressures, these  curves  help  narrow  the  possible 
range  of  radiant  flux  densities  and  will  be  more  useful  if  and  when 
more  becomes  known  about  temperature  profiles  at  high  pressure. 

The  sensitivity  of  radiant  heat  flux  to  temperature  was  inves¬ 
tigated  for  a  sequence  of  profiles  of  the  form  of  Eq.  3.  3  (with  x=3) 
in  which  centerline  temperature  was  varied  from  12,  000  to  15,  000  K, 
Fig.  4.3.  Over  the  range  calculated  the  radiant  heat  flux  increases 
approximately  as  the  8.4  power  of  arc  column  centerline  temperature 
at  100  atm  and  below, 

4.4  SENSITIVITY  TO  CHANGES  IN  OUTER  RADIUS 

Weber's  estimates  of  constricted  arc  column  growth  (Ref.  12) 
indicate  that  the  arc  may  not  entirely  fill  the  constrictor  of  the  arc 
heater  used  in  the  experimental  work.  Temperature  could  not  be 
measured  over  the  entire  column,  so  for  one  of  the  measured  profiles, 
Fig.  5.4,  two  different  outer  edge  parabolas  were  matched  to  the 
polynomial  which  was  fitted  to  the  data  (in  this  case  also  a  parabola). 
The  part  of  the  profile  for  which  data  existed  was  unchanged  and  the 
sensitivity  of  the  calculation  to  changes  in  the  outer  part  of  the  pro¬ 
file  was  checked  in  this  way.  Comparison  made  on  the  basis  of  heat 
flux  at  the  constrictor  wall  (0.5  cm  radius)  indicated  that  exact  arc 
column  diameter  made  little  difference.  Table  4.2. 

Because  only  the  outer  few  annuli  were  changed  this  result 
does  not  contradict  Fig.  4.2  where  the  entire  profile  was  systemat¬ 
ically  varied. 
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tcl  black  body  limit 


RADIANT 


-  LF  +  HF 

- LF  see  sect.  4.1 

PROFILE  EQ  3.3 
T  ref.  «  2000°  K  X-3, 
8  =  0,5  cm 


12000  13000  14000  15000 

ARC  CENTERLINE  TEMPERATURE  (°K) 

Figure  4,  3  Temperature  Level  Dependence  of  Arc  Radiant  Heat  Flux, 
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Pressure 


Radiant  Flux  Density  at  0.5  cm  radius  (kw/cm  ) 


(atm) 


6  =  .  35  cm 

6  =  .  5  cm 

HF 

.44 

.62 

LF 

.  18 

.20 

SUM 

.62 

.82 

HF 

.88 

.93 

LF 

.55 

.59 

SUM 

.  1.43 

1.  52 

HF 

1.87 

1.65 

LF 

4.62 

4.88 

SUM 

6.49 

6.53 

HF 

2.5 

3.2 

LF 

45.  1 

45.6 

SUM 

47.6 

48.8 

Table  4.  2  -  Radiant  Flux  Density  at  Constrictor  Boundary 
2.4  atm  455  amp  measured  profile 


4.5  QUASI-THIN  MODEL  AND  FREQUENCY  DISTRIBUTION 

The  quasi-thin  model  was  obtained  by  setting  all  transmis¬ 
sivities  equal  to  unity  in  Eq.  2.26.  Because  the  computation  of 
emissivities  takes  into  account  self  absorption  within  each  annulus, 
the  quasi-thin  model  will  yield  a  radiant  heat  flux  known  to  be  closer 
to  the  true  value  than  the  optically  thin  model.  The  extent  of  the 
difference  between  quasi-thin  and  optically  thin  depends  on  the  mag¬ 
nitude  of  the  dimensionless  absorption  lengths  associated  with  each 
annulus,  so  choice  of  number  of  annuli  will  have  an  effect  on  the 
quasi-thin  result.  The  quasi-thin  result,  however,  is  a  very  simple 
extra  step  for  the  computer  and  provides  a  useful  indication  of  the 
importance  of  self-absorption. 

As  mentioned  in  section  3.  1,  an  option  in  the  program  permits 
automatic  plotting  of  radiant  flux  density  in  each  wavenumber  interval, 
both  quasi-thin  and  self- absorbing.  Fig.  4.4.  In  general  the  higher 
the  intensity  of  radiation  the  greater  the  difference  between  self¬ 
absorbing  and  quasi-thin. 
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Figure  4.4  Frequency  Distribution  of  Radiant  Heat  Flux.  Arc  Profile  Measured 
at  2.  4  atm,  455  amp.  Computed  at  100  atm. 
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4.6  TWO  TEMPERATURE  MODEL 

In  order  to  gain  some  idea  of  the  effect  of  a  body  of  heated 
plasma  in  reducing  the  intense  radiation  from  the  arc  column  proper, 
a  two  temperature  cylindrical  model,  which  is  compatible  with  that 
discussed  in  section  2,  was  investigated.  Figure  4.5  presents 
results  with  the  arc  at  16,000°K  and  the  surronding  plasma  at  6,  000°K 
with  the  dimensions  investigated  shown  in  the  inset  sketch.  As  pressure 
is  increased,  the  shielding  effect  of  the  lower  temperature  surrounding 
plasma  may  be  sufficient  to  cause  the  total  radiant  flux  to  pass  through 
a  maximum  and  then  decrease  somewhat.  As  pressure  is  increased 
the  energy  radiated  by  the  cooler  surrounding  gas  to  the  walls  will 
continue  to  increase  and  the  large  equivalent  radiating  surface  tends 
to  magnify  the  increase.  However,  the  fraction  of  energy  reabsorbed 
from  the  arc  column  will  also  increase  and  the  possibility  that  the 
latter  effect  could  result  in  a  net  decrease  in  radiant  flux  to  the  walls 
in  confirmed  by  Fig.  4.5.  The  calculation  waB  repeated  with  the  arc 
at  15,  000°K  in  order  to  include  the  high  frequency  radiation.  Fig.  4.6. 

It  must  be  emphasized  that  the  existence  of  Buch  an  effect  is 
controlled  by  geometry,  operating  conditions  and  arc  behavior  at  high 
pressures  --  all  of  which  affect  the  temperature  distribution  and  its 
{here  neglected)  pressure  dependence. 
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Figure  4,  5  Reduction  in  Intensity  of  Radiation  from  a  High  Pressure  by  Surrounding 
Plasma  -  Two  Temperature  Model  with  Arc  at  16,  000  K. 
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s  s  i 


RADIATION  FROM  ARC 
COLUMN  REACHING 
THE  WALL 


100  200 


500  1000 


Figure  4,6  Two  Temperature  Model  Including  High  Frequency  Radiation 
with  Arc  at  15,  000°K. 


ll'£9'aX'D03V 


AEDC-TR-65.11 


SECTION  V 

EXPERIMENTAL  STUDY 


5.1  OBJECTIVES 

The  objectives  of  the  experimental  work  were: 

1)  measurement  of  temperature  profiles  for  use  as  an 
input  to  the  analysis. 

2)  measurement  of  total  radiant  heat  flux  to  compare  with 
results  of  the  analysis. 

5.2  ARC  HEATER 

The  Gerdien  type  arc  heater  used  for  this  study  is  shown  in 
Fig.  5.  1.  This  heater  produces  a  uniform  uncantaminated  air  plasma 
in  a  test  section  (constricted  region)  instrumented  with  voltage  probes, 
a  spectro graphic  observation  window  and  a  shuttered  channel  for 
radiation  measurements.  The  heater  has  been  more  completely  des¬ 
cribed  in  Ref.  8.  The  constrictor  section  can  be  replaced  by  a  some¬ 
what  shorter  non- constricting  section  which  has  a  window  for  photo¬ 
graphic  observation. 

Available  power  supply  voltage  limits  maximum  operating 
pressure  in  the  unconstricted  mode  to  15  atmospheres  although  the 
heater  is  capable  of  withstanding  higher  pressures.  The  greater 
length  of  the  constrictor  section  reduced  the  voltage  limitation  on 
pressure  to  about  5  atmospheres. 

5.3  TEMPERATURE  PROFILES  AND  CENTERLINE  TEMPERATURES 

Two  temperature  profiles  were  obtained  at  one  atmosphere 
and  one  at  2.4  atmospheres.  These  are  shown  in  Figs.  5.2,  5.3, 
and  5.4.  The  reason  for  two  fitted  curves  in  Fig.  5.4  was  discussed 
in  Section  4. 4. 

Relative  intensities  of  several  emission  lines  were  used  with 
Abel's  integral  equation  to  establish  arc  column  temperature  gradients. 
A  description  of  the  technique  is  included  as  Appendix  E.  Additional 
measurements  of  arc  centerline  temperature  were  made  at  2.2 
atmospheres  and  several  current  levels  as  Bhown  in  Fig.  5.  5.  These 
are  based  on  average  intensities  in  a  narrow  slice  through  the  center 
of  the  arc.  Therefore  evaluation  cannot  make  use  of  Abel's  integral 
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J _ 1 _ I _ I _ I _ I 

.1  .2  .3  .4  .5  .6 

r  (cm) 

Spectrographically  Measured  Profile,  1  atm,  585  amp. 


45 


EMPERATURE  (°K) 


A  EDC-TR-65-1 1 


16000 


ANALYTIC  APPROXIMATION 


14000 


12000 


10000 


8000 


6000 


2.4  atm 
455  amp 

O  RATIO  OF 
X7774  and  X7947 


ARC  COLUMN  RADIUS  (cm) 


Figure  5.  4  Spectrographically  Measured  Profile,  2.  4  atm,  455  amp 
with  Curve  Fit  to  Two  Different  Boundary  Conditions. 
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to  allow  for  the  cooler  gas  ahead  of  and  behind  the  centerline  region 
and  the  resultant  temperatures  are  slightly  lower  that  the  "true" 
centerline  temperature  as  obtained  by  measurements  over  a  signi¬ 
ficant  fraction  of  the  profile. 

Spectrographic  data  were  not  obtained  at  higher  pressures 
because  the  measurements  depend  on  observed  spectral  lines  being 
nearly  optically  thin.  The  line  reversal  technique  which  is  applicable 
at  high  pressures  requires  of  the  order  of  70%  absorption  and  is  there¬ 
fore  most  applicable  at  pressures  substantially  above  5  atm. 

SA  RADIANT  HEAT  FLUX  MEASUREMENTS 

Three  types  of  gages  were  considered  for  measuring  radiant 
heat  flux  while  avoiding  convective  heating.  The  first  gage  was  a 
cavity  type,  recessed  in  a  channel  5  cm  from  the  arc  column.  The 
channel  eliminated  convection  and  conduction  without  the  necessity 
of  a  transparent  window,  but  temperature  rise  of  the  cavity  was 
barely  detectable.  At  the  other  extreme  a  thermopile  proved  too 
sensitive. 

The  gage  actually  used.  Fig.  5.6,  was  a  cylindrical  copper 
slug  mounted  with  one  base  coated  with  lampblack  and  covered  by 
a  l/16Mthick  quartz  window.  Remaining  surfaces  were  thermally 
insulated  by  a  teflon  sleeve  which  also  held- the  window  in  place. 
Temperature  rise  in  the  slug  was  measured  by  a  chromel-alumel 
thermocouple  directly  connected  to  a  recording  oscillograph. 

A  shutter  shielded  the  gage  until  the  arc  had  stabilized  at  the 
desired  conditions.  Then  the  shutter  was  opened  and  the  gage  moved 
rapidly  forward  to  one  of  several  preselected  positions  in  the  channel. 
In  the  first  position  the  gage  surface  was  recessed  1/2"  from  the 
constrictor  wall.  Fully  extended  the  gage  surface  was  recessed  1/16" 
so  that  the  exposed  surface  of  the  quartz  was  tangent  to  the  constrictor 
wall. 


Radiant  heat  flux  to  the  gage  was  calculated  from  the  rate  of 
temperature  rise  of  the  copper  slug,  assuming  constant  specific  heat 
and  no  heat  losses.  A  typical  data  trace  is  shown  in  Fig,  5.7. 
Radiation  ot  the  constrictor  wall  can  be  inferred  by  plotting  heat 
flux  data  at  several  points  along  the  channel  and  extrapolating  to  zero 
distance. 
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to  allow  for  the  cooler  gas  ahead  of  and  behind  the  centerline  region 
and  the  resultant  temperatures  are  slightly  lower  that  the  "true" 
centerline  temperature  as  obtained  by  measurements  over  a  signi¬ 
ficant  fraction  of  the  profile. 

Spectrographic  data  were  not  obtained  at  higher  pressures 
because  the  measurements  depend  on  observed  spectral  lines  being 
nearly  optically  thin.  The  line  reversal  technique  which  is  applicable 
at  high  pressures  requires  of  the  order  of  70%  absorption  and  is  there¬ 
fore  most  applicable  at  pressures  substantially  above  5  atm. 

5.4  RADIANT  HEAT  FLUX  MEASUREMENTS 

Three  types  of  gages  were  considered  for  measuring  radiant 
heat  flux  while  avoiding  convective  heating.  The  first  gage  was  a 
cavity  type,  recessed  in  a  channel  5  cm  from  the  arc  column.  The 
channel  eliminated  convection  and  conduction  without  the  necessity 
of  a  transparent  window,  but  temperature  rise  of  the  cavity  was 
barely  detectable.  At  the  other  extreme  a  thermopile  proved  too 
sensitive. 

The  gage  actually  used,  Fig.  5,  6,  was  a  cylindrical  copper 
slug  mounted  with  one  base  coated  with  lampblack  and  covered  by 
a  l/16"thick  quartz  window.  Remaining  surfaces  were  thermally 
insulated  by  a  teflon  sleeve  which  also  held  the  window  in  place. 
Temperature  rise  in  the  slug  was  measured  by  a  chromel-alumel 
thermocouple  directly  connected  to  a  recording  oscillograph. 

A  shutter  shielded  the  gage  until  the  arc  had  stabilized  at  the 
desired  conditions.  Then  the  shutter  was  opened  and  the  gage  moved 
rapidly  forward  to  one  of  several  preselected  positions  in  the  channel. 
In  the  first  position  the  gage  surface  was  recessed  1/2"  from  the 
constrictor  wall.  Fully  extended  the  gage  surface  was  recessed  1/16" 
so  that  the  exposed  surface  of  the  quartz  was  tangent  to  the  constrictor 
wall. 


Radiant  heat  flux  to  the  gage  was  calculated  from  the  rate  of 
temperature  rise  of  the  copper  slug,  assuming  constant  specific  heat 
and  no  heat  losses.  A  typical  data  trace  is  shown  in  Fig.  5.7. 
Radiation  ot  the  constrictor  wall  can  be  inferred  by  plotting  heat 
flux 'data  at  several  points  along  the  channel  and  extrapolating  to  zero 
distance. 
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Allowances  must  be  made  for  losses  through  the  window  and 
at  the  gage  surface*  The  percentage  of  light  reflected  by  a  glass 
surface  of  refractive  index  1.55  (about  that  of  quartz)  is  tabulated 
as  a  function  of  angle  of  incidence  in  Ref.  13  and  a  graphical  inte¬ 
gration  based  on  this  tabulation  indicates  that  approximately  92%  of 
the  incident  energy  is  transmitted  when  the  gage  is  fully  in.  Taking 
the  absorptivity  of  lampblack  as  .95,  the  radiant  energy  measured  by 
the  gage  at  the  arc  chamber  wall  must  be  increased  by  a  factor  of 
1,  14.  This  correction  is  approximate  but  extrapolation  to  the  arc 
column  wall  leads  to  considerable  uncertainty  in  any  event.  The 
transmission  characteristic  of  type  101  fused  quartz.  Fig.  5.8,  coin¬ 
cides  very  nearly  with  the  Low  Frequency  Range  computer  output  so 
comparison  will  be  made  with  LF  curves. 

An  additional  uncertainty  introduced  by  the  quartz  window  is 
the  possible  change  in  transmission  characteristic  when  the  glass  is 
exposed  to  the  high  temperature  arc.  The  transmission  characteristic 
in  the  visible  range  was  checked  for  a  few  windows  using  a  standard 
lamp  and  photomultipliers.  The  ratio  of  photomultiplier  outputs  with 
and  without  the  window  inserted  gave  the  fraction  of  energy  transmitted 
(i.  e. ,  energy  neither  reflected  nor  absorbed).  A  clean  window  gave 
values  from  .  89  to  .  90,  Transmission  of  windows  which  has  been  ex¬ 
posed  to  the  arc  ranged  from  .40  to  ,  79.  The  time  at  which  the  changes 
took  place  is  unknown  so  no  correction  has  been  made  for  this  effect. 

Results  are  shown  in  Table  5.  1  and  Fig.  5.  9.  For  the  first 
series  the  insulator  used  to  jacket  the  copper  slug  was  melamine, 
and  the  quartz  window  was  fastened  to  the  melamine  with  an  adhesive. 
This  was  not  entirely  satisfactory  since  the  glass  did  not  always  re¬ 
main  in  place.  In  the  second  series  tabulated  the  teflon  sleeve  arrange¬ 
ment  shown  in  Fig.  5.  6  proved  to  be  much  more  reliable. 

The  effect  of  the  quartz  window  in  shielding  the  gage  from  con¬ 
duction  is  evident  in  the  comparison  of  table  5.2. 

Figure  5.  10  is  a  comparison  of  calculations  based  on  two 
measured  profiles  and  the  results  of  direct  measurement.  Complete 
confirmation  of  the  calculation  model  and  absorption  coefficient  data 
requires  measurements  at  substantially  higher  pressures  but  at  least 
the  magnitude  and  trend  are  consistent. 
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Figure  5.8  Average  Transmission  of  General  Electric  Type  101  Quartz  Excluding 
Surface  Reflection  Losses  (Th!ckness  =  l/l6M). 
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Figure  5.  10  Comparison  of  Computations  Using  2  Measured  Profiles 
with  Direct  Radiant  Flux  Density  Measurements 


54 


AEDC-TR-65-1 1 


Distance  of 
Gage  Face 
from  Arc 
Column 
(in) 

Arc 

Current 

(amp) 

Arc 

Voltage 

(volts) 

Rate  of 
Temperature 
i  Rise 

(°F/  sec) 

Mass 

of 

Slug 

(gm) 

Radiant 
Flux 
Density 
(w/ cm^) 

First  series,  P=2.  2  atm 

1.  25 

460 

350 

6.  86 

1.  64 

7.7 

1.25 

400 

410 

5.  8 

l.  64 

6.5 

1.  25 

430 

410 

6.  7 

1.  64 

7.  5 

1.  25 

405 

390 

5.8 

1.  64 

6.5 

1.  25 

430 

415 

8.  1 

1.  64 

9.1 

1.  25 

525 

355 

6.  9 

1.75 

8.  2 

.  50 

450 

390 

11.  86 

1.  64 

13.  3 

.50 

440 

350 

21. 0 

1.  64 

23.  6 

.  50 

425 

405 

11.  60 

1.  64 

13.  0 

.28 

440 

360 

50.  2 

1.  64 

56.4 

.  22 

500 

365 

76.4 

1.  64 

85.  7 

.13 

500 

356 

172.  0 

1.  53 

179.5 

.13 

550 

340 

175.  0 

1.53 

183.  0 

.13 

545 

345 

94.5 

1.75 

112. 

.13 

555 

340 

109.0 

1.75 

129. 

.13 

555 

330 

99.0 

1.  75 

117.0 

.06 

485 

372 

179.0 

1.  72 

210.  5 

.  06 

550 

350 

136.  0 

1.  75 

161.  0 

Second  series,  F=2,  2  atm 

.  06 

400 

430 

249 

1.  67 

284 

.  06 

480 

410 

244 

1.  67 

278 

.  06 

440 

408 

252 

1.  67 

292 

.  06 

450 

410 

266 

1.  67 

302 

.13 

450 

405 

156 

1.  67 

178 

.13 

455 

415 

185 

1.  67 

211 

.25 

425 

415 

123 

1.  67 

141 

Second  series,  P=5.  7  atm 

.  06 

600 

505 

479 

1. 76 

575 

.  06 

610 

480 

385 

1.  76 

454 

Table  5.1-  Radiant  Heat  Flux  to  Slug  Gage 
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Distance  of 
Gage  Face  from 
Arc  Column 
(cm) _ 


.  559 
.159 
.  159 


Table  5.  2 


Radiant  Flux 
Density 
(watts /cm^) 


Total 

Flux  Density 
(No  Window  on  Gage) 
(watts  /cm^) _ 


85.  7 


214.  0 


210.5 


636,  0 


161.  0 


648.  0 


Comparison  of  Total  and  Radiant  Heat  Flux 
at  a  Pressure  of  Two  Atmospheres 
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SECTION  VI 

APPLICATION  TO  ARC  HEATER  CONFIGURATIONS 


6.1  GENERAL 

A  variety  of  configurations  for  heating  plasma  by  means  of  an 
electric  arc  are  in  use  today.  Two  of  the  more  common  basic  arc 
heater  geometries  are  the  constricted  arc  and  the  ring  electrode 
rotating  arc.  Application  of  the  model  developed  in  Section  2  applied 
fairly  well  to  the  constricted  arc  since  axial  symmetry  is  maintained 
in  most  of  the  heater  and  some  measurements  of  arc  column  tem¬ 
perature  profiles  are  available. 

Analysis  of  the  ring  electrode  rotating  arc  is  considerably 
more  difficult  because  arc  temperature  measurements  are  not  avail¬ 
able  and  arc  behavior  is  influenced  by  such  new  factorB  as  the  ratio  of 
ring  diameter  to  separation  distance,  the  magnitude  and  direction  of 
self -induced  or  externally  applied  magnetic  driving  fields,  and  the 
magnitude  and  direction  of  cold  air  flow  transverse  to  the  arc  column. 

Under  some  circumstance s  for  instance,  the  rotating  arc  may 
bow  inward  to  the  extent  that,  except  near  the  electrodes,  it  is 
stationary  on  the  axiB  of  symmetry,  Fig.  6.  1.  Fastax  films  taken  of 
a  rotating  arc  between  ring  electrodes  separated  by  a  distance  of 
about  4  ring  radii  tend  to  bear  out  this  possibility.  The  column  was 
photographed  through  a  3/8"  wide  x  1/8”  high  port  midway  between 
electrodes.  At  pressures  above  approximately  10  atmospheres  the 
the  column  was  visible  and  stationary,  centered  in  the  port.  At  lower 
pressures,  down  to  about  2  atmospheres,  changes  in  brightness 
indicated  that  the  column  was  moving  past  the  port.  At  2  atmospheres 
pressure  light  from  the  port  was  uniform.  These  observations  fit  the 
hypothesis  that  at  the  highest  pressure  only  the  parts  of  the  arc  near 
the  electrodes  are  moving,  while  the  body  of  the  column  is  aligned 
with  the  axis,  while  at  low  pressure  the  entire  column  may  be  moving 
too  fast  for  the  Fastax  camera  to  "stop"  the  motion. 

If  the  magnetic  field  which  causes  the  arc  to  rotate  is 
sufficiently  strong,  of  the  order  of  20  kilogauss,  the  arc  can  be 
magnetically  diffused  (Ref.  14)  thus  radically  altering  its  geometry 
and  radiation  characteristics.  A  transverse  flow  created  by  a  rapidly 
moving  arc  column  may  destroy  the  axial  symmetry  of  the  column, 
probably  tending  to  create  a  turbulent  If  thiB  tail  merge b  with  the 
leading  "edge"  of  the  column  a  diffuse  operating  mode  also  results. 
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RING  ELECTRODE 
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HIGH  PRESSURE 
CENTER  OF  ARC 
STEADY. 


Figure  6.  1  Ring  Electrode  Rotating  Arc  With  Separation  Distance  of 
4  Ring  Radii. 
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6.2  SPECIFIC  CONFIGURATION 

An'outline  sketch  of  the  arc  heater  of  immediate  interest  is 
shown  in  Fig.  6.2.  It  is  designed  to  operate  at  high  pressure  on  3 
phase,  60  cps,  alternating  current  with  the  arcs  kept  in  rapid  motion 
by  a  steady,  externally  applied,  magnetic  field.  Cold  air  flows 
radially  inward  through  the  arcs,  then  axially  to  the  nozzle. 

Magnetic  field  strength  is  5  kilogauss,  which  is  probably  not 
sufficient  to  cause  magnetic  diffusion,  so  the  arc  will  be  treated  as  a 
cylindrical  column.  Arc  length  cannot  easily  be  established  because, 
while  electrode  separation  is  only  about  one  fourth  that  in  the  observa¬ 
tions  described  in  Section  6.  1,  pressure  can  be  higher  and  there  is  a 
radially  inward  velocity  component  which  will  probably  cause  the  arc 
to  bow  inward.  Some  estimated  operating  conditions  at  two  pressures 
are  shown  in  Table  6.  1. 

TABLE  6.  1.  AEDC  20  MEGAWATT  ARC  HEATER  SPECIFICATIONS 

Dimensions 


Chamber  Diameter  (in. )  4 

Insulator  Diameter  (in. )  6 

Electrode  Spacing  (in.  ) 

center  to  center  4.  5 

Minimum  Gap  1 . 5 

Total  Chamber  Length  (in.  }  17.  5 


Low 

Pressure 


Pressure  (psia)**  75 

(atm)  5.  1 

Arc  Voltage  (volts  zero  to  mean)**  1000 

Arc  Current  (amperes  rms/phase)**  4000 

Power  Input  (MW)**  11 

Exit  Enthalpy  (Btu/lbm)**  6100 

(h/RTo)  180 

Exit  Temperature  (°K)  6300 

Mass  Flow  (lbm/Bec)**  1.20 

Estimated  Rotational  Speed  (rev/sec)**  1400 


♦Scaled  from  Westinghouse  Drawing  AC09084 
**  Westinghouse  letter  of  1 1  September  1964 


High 

Pressure 

3750 

255 

2000 

3500 

19 

7200 

212 

8000 

1.0 

1300 
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6.3  ARC  COLUMN  MODEL 

Arc  size  and  temperature  distribution  are  not  known,  so  some 
admittedly  crude  assumptions  were  made  in  order  to  arrive  at  an 
estimate.  The  arc  was  assumed  to  be  cylindrical  at  a  uniform  tem¬ 
perature  with  air  flow  through  it  limited  by  aerodynamic  choking  or 
by  the  total  flow,  whichever  was  least.  The  choked  flow  function  for 
flow  through  an  arc  has  been  computed  by  Weber  (Ref.  15)  and  is 
shown  in  Fig.  6.  3.  For  a  known  or  assumed  voltage  gradient  and 
current  there  is  then  a  unique  combination  of  arc  radius  and  tem¬ 
perature  which  satisfies  both  an  energy  balance  and  Ohm's  Law.  In 
view  of  the  approximate  nature  of  the  data  and  model,  electrode  drops 
{of  the  order  of  tens  of  volts)  and  back  EMF  due  to  interaction  with  the 
externally  applied  magnetic  field  {a  few  volts)  were  neglected  in  the 
estimates  of  voltage  gradient.  A  further  implicit  assumption  is  that 
an  AC  arc  may  be  treated  by  DC  methods  and  that  the  RMS  value  of 
current  gives  acceptable  mean  values.  In  brief,  both  space  and  time 
variation  of  arc  temperature  (and  bulk  gas  temperature  too)  were 
ignored. 

6.4  LOW  PRESSURE  OPERATION 


The  gas  is  nearly  optically  thin  at  5  atmospheres  so  that  total 
radiance,  J,  as  computed  by  Breene  and  Nardone  (Ref.  2)  was  used  to 
computed  radiant  heat  fluxes.  The  energy  balance  for  a  unit  length  of 
arc  column  was  written 


<r  fffl2 


e 


p  (2d)  (h  -  h  )  -  (4tt  J  )  (tt  6  )  =  0 
a  o  a 


(6.1) 


where  p  is  pressure,  h  is  the  enthalpy  of  the  incoming  air,  is 
enthalpy  at  arc  temperature,  6  is  arc  radius  and  (— ^)  is  the  choked 
flow  function.  Fig.  6.3.  Peng  and  Pendroh  (Ref.  1 6)  ^were  the  source 
for  electrical  conductivity  Ce.  Equation  6.  1  can  be  considered  a 
function  of  radius  and  arc  temperature,  Fig.  6.4,  since,  at  a  given 
pressure,  C7e,  J,  ( )  and  ha  are  all  functions  of  temperature. 

Ohm's  law  in  the  form 

E  =  - - L_-  (6.2) 

cr  it  6 

e 

is  also  shown  pn  Fig.  6.4  for  several  values  of  E.  The  values  selected 
are  reasonable  at  the  low  pressure.  Experience  with  rotating  arcs 
indicates  that  the  column  radius  should  be  less  than  about  0.  5  cm  and 
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that  60  v/cm  is  a  reasonable  maximum  voltage  gradient.  This  is  a 
local  column  gradient  which  cannot  be  obtained  by  simply  dividing  the 
voltage  between  electrodes  by  electrode  spacing.  Turbulence  locally 
twists  and  turns  the  column  extending  its  actual  length  to  several 
times  the  distance  between  electrodes.  The  arc  length  used  for 
calculations  at  5  atm  is  then  obtained  by  dividing  the  applied  voltage 
by  the  voltage  gradient. 

Bulk  enthalpy,  h^,  of  plasma  leaving  the  arc  chamber,  was 
estimated  from  another  energy  balance.  Consider  3  arcs,  each  of 
length  L,a,  in  a  cylindrical  chamber  of  length  Lc  and  radius  Rc  and 
assume  that  the  only  loss  is  the  bulk  gas  radiation  at  bulk  (exit) 
enthalpy.  Then 

p(26)(3L)h  -  Wh  =  J.  (ff R  2  L  )  (6.3) 

3.  b  D  C  C 

where  subscript  b  refers  to  bulk  gas  and  La  is  the  length  of  a  single 
arc.  The  choked  flow  function  and  the  total  radiance  are  known 
functions  of  bulk  enthalpy  (or,  equivalently,  bulk  temperature)  so 
Eq.  6.  3  uniquely  determines  h^.  Results  are  presented  in  Section 
6.  6  along  with  those  for  the  high  pressure  case. 

6.5  HIGH  PRESSURE  OPERATION 

Essentially  the  same  analysis  was  applied  to  the  255  atmos¬ 
phere  case  except  that  the  gas  can  no  longer  be  considered  optically 
thin  and  program  ARCRAD  was  used  to  calculate  radiation.  For  the 
conditions  given,  the  entire  flow  can  pass  through  the  arc  column 
without  choking,  so  the  arc  column  energy  balance  analogous  to 
Eq.  6.  I  was  written 

2  W(h  -h) 

- 1 - - (q/A)  2  IT  6  ~  0  (6.4) 

<T  v  fl  3  a  a 

e 

where  (q/A)a  is  the  radiant  flux  density  at  the  arc  column  boundary. 
Fig.  6.5. 

Again  a  unique  combination  of  a  radius  and  arc  temperature 
exist  which  simultaneously  satisfy  Eq.  6.4  and  Ohm's  law  (Eq.  6.2) 
for  an  established  value  of  E. 

Raezer,  et,  al.  ,  Ref.  (17)  indicate  that  voltage  gradient  varies 
as  the  0.28  power  of  the  pressure.  Applying  this  to  the  60  v/cm 
gradient  at  5  atmospheres  gives  174  v/cm  and  an  arc  length  of  4.  5 
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inches.  The  resulting  arc  radius  is  quite  large,  probably  much  too 
large,  and  the  arc  temperature  and  enthalpy  are  therefore  low  since 
the  general  trend  of  the  curve  is  similar  to  Fig,  6.4.  For  this  reason 
results  are  also  included  for  an  arbitrarily  assumed  arc  length  of 
2  inches. 


Exit  enthalpy  was  again  estimated  from  a  simplified  arc 
chamber  energy  balance.  The  bulk  gas  was  taken  to  be  a  cylinder  of 
arc  chamber  dimensions  at  a  uniform  temperature  corresponding  to 
exit  enthalpy,  with  a  smaller  cylinder,  the  arc  column,  moving  rapidly 
around  its  circumference.  Both  the  bulk  gas  and  the  arc  column  radiate 
to  the  arc  chamber  wall  but  a  part  of  the  arc  radiation  is  now  re¬ 
absorbed  by  the  bulk  gas. 


In  cross  section  this  geometry  corresponds  very  closely  to 
radiation  of  a  surface  element  into  a  cylinder,  so  the  Nusselt  solution 
for  absorptivity  of  a  cylinder  can  be  applied.  By  Lambert’s  cosine  law 
the  equivalent  arc  surface  area  per  unit  length  radiating  inward  is  equal 
to  the  arc  column  diameter  so  that  qab»  the  power  reabsorbed  by  the 
bulk  gas,  is  given  by 


qab 


€  V  a 
a,n  a,n  b,n 


(6.5) 


Program  ARCRAD  was  adapted  to  the  computation  of  q  . 

ab 

Neglecting  all  losses  but  radiation  and  assuming,  optimistically, 
that  radiation  from  the  bulk  gas  takes  place  at  bulk  exit  enthalpy,  the 
equivalent  of  Eq.  6.  3  is 


W(h  -  h  } 
a  b 


+  qab  *  qt 


(6.6) 


where  power  radiated  by  the  bulk  gas, 


qbis 


given  by 


qb  =  (q/A)b  (2ffRc)  (Rc  +  L^  (6.7) 

and  {q/A)^  was  computed  using  program  ARCRAD.  Strictly  speaking, 
absorption  lengths  for  radiation  from  the  ends  of  the  cylinder  are 
slightly  different  than  for  the  cylindrical  surface,  but  the  difference 
is  of  the  order  of  5%  and  the  ends  are  only  a  small  fraction  of  the 
total  surface  anyway  so  no  distinction  is  made.  It  was  convenient  to 
solve  for  T^  graphically,  Fig,  6.  6,  plotting  Eqs.  6.  6  and  6,  7  as 
functions  of  temperature  because  both  and  h^  are  non-analytic. 
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POWER  TO  BULK  GAS 
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l _ I _ I _ I 

5000  6000  7000  8000 

BULK  TEMPERATURE  (®K) 


6  Arc  Chamber  Energy  Balance  at  255  atmospheres  with 
an  Input  power  of  19  megawatts. 
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Two  curves  are  shown  for  q^,  one  assuming  that  the  chamber 
walls  are  black  and  the  other  applying  the  correction  factor  discussed 
in  Section  2.  7.  The  arc  itself  is  small  compared  to  the  chamber  so 
no  correction  is  needed,  but  reabsorbtion  by  the  bulk  gas  can  be 
significant.  In  evaluating  the  correction  factor  for  the  bulk  gas  a  g 
was  taken  as  the  ratio  of  (q/A)^  to  black  body  radiation,  both  at 
7500  °K, 

If  Pioss  is  defined  as  that  part  of  the  input  power  which  is 
absorbed  by  processes  other  than  net  heating  of  the  incoming  air  and 
Pin  is  the  total  input  power  then 


P  £  P 

loss  in 


W( h  -  h  ) 
b  o 


(6.8) 


A  curve  representing  the  equality  is  also  plotted  in  Fig.  6.  6  and  its 
intersection  with  the  radiation  loss  curve  places  an  upper  bound  on 
bulk  temperature. 

6.6  RESULTS  AND  DISCUSSION 

Results  of  calculations  are  summarized  below. 
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TABLE  6.2.  AEDC  20  MEGAWATT  ARC  HEATER,  ESTIMATED  ARC 
CONDITIONS  AND  RADIANT  FLUXES 


Pressure  (atm) 

Voltage  gradient  (V/cm) 

5 

60 

255 

390 

255 

174 

Arc  Length  (cm) 

15.  7 

5. 

1 

11.5 

(in) 

6.2 

2. 

0 

4.  5 

Arc  Radius  (cm) 

.48 

9 

45 

.  78 

Arc  Temperature  (K) 

16,  000 

11, 

000 

10,500 

Wall  Emissivity 

.57 

1.0 

.  57 

1.0 

Exit  Temperature  (K) 

6,  300 

7500 

7200 

7400 

7000 

Exit  Enthalpy  (Btu/lbm) 

6,100 

6300 

5900 

6200 

5600 

Arc  Radiation  (mw) 

Total  from  arc 

3.6 

.9 

.9 

3.  6 

3.6 

To  Bulk  gas* 

0 

.2 

.2 

.  6 

.6 

To  Wall* 

3.6 

.  7 

.7 

3.0 

3.  0 

Bulk  gas  radiation  (mw) 

.05 

11.  0 

11.7 

9.7 

10.  3 

Total  Radiation  to  Wall 

3.65 

11.7 

12.4 

12.7 

13.  3 

Average  Radiant  Flux 
Density**  (kw/cm^) 

2.3 

7.4 

7.8 

8.0 

B.  4 

*  at  7500°K 

**  Based  on  a  cylinder  17.  5  inches  long  by  4  inches  diameter. 

At  5  atmospheres  the  bulk  gas  radiation  loss  is  relatively 
small  and  the  decrease  in  enthalpy  from  arc  to  bulk  gas  is  primarily 
the  result  of  mixing  with  cold  gas  which  did  not  pass  through  the  arc. 

At  255  atmospheres  bulk  gas  radiation  alonge  can  consume 
practically  all  of  the  power  input.  Results  presented  in  Fig.  6.  6  and 
Table  6.  2  tend  to  be  optimistic  in  that  conduction  losses  have  been 
neglected  and  bulk  gas  radiation  has  been  assumed  to  take  place  at 
exit  temperature  rather  than  somewhere  between  arc  temperature  and 
bulk  temperature. 
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Total  radiation  to  the  wall  includes  radiation  directly  from  the 
arcs.  The  arcs  would  tend  to  create  local  hot  spots  except  that  they 
are  in  rapid  motion  --  making  about  10  revolutions  for  each  half  cycle 
of  input  voltage  --so  the  heat  flux  per  unit  length  from  the  arc  was 
considered  as  uniformly  distributed  around  the  arc  chamber  periphery. 
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SECTION  VII 
CONCLUSIONS 

At  high  pressures  there  is  no  simple  way  to  predict  the  radiant 
heat  transfer  from  an  arc  column.  The  computational  model  des¬ 
cribed  herein  provides  very  tight  upper  and  lower  limits  if  the  basic 
geometry  assumptions  of  the  model  are  satisfied  and  if  the  temper¬ 
ature  distribution  is  known.  Reasonable  limiting  assumptions  about 
temperature  profiles  can  be  useful  in  bracketing  the  desired  answer 
but  one  must  know  the  temperature  level  either  by  measurement  or 
calculation. 

Incorporation  of  radiation  including  self -absorption  into  a 
solution  of  the  energy  equation  such  as  Weber  has  done  {Ref.  11)  for 
the  optically  thin  case,  is  the  next  logical  step  in  the  study  of  this 
problem.  The  temperature  profile  would  then  be  a  part  of  the  result 
rather  than  a  required  input.  A  coordinated  experimental  program  to 
measure  radiant  heat  flux  and/or  temperature  profiles  at  high  pressure 
would  also  be  very  desirable. 
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APPENDIX  A 

CALCULATION  OF  ENERGY  FRACTION 

For  a  given  wavenumber  interval,  An  =  n.  -n  ,  let  1)  be  the 
fraction  of  the  total  energy  flux  radiated  by  a  black  body  at  a  given 
temperature  which  falls  within  this  interval 


tj=D 


1 


where 


CD 


o 


(A.  1) 


(A.  2) 


The  function  D  is  available  in  tabular  form,  (Ref.  18)  but,  for  computer 
application,  direct  generation  of  T)  by  series  expansion  was  desired. 
Two  series  were  used. 


1.  Pivovonsky  and  Nagel  (Ref.  18)  gives  a  series  expansion 
for  D  as 


D  =  1 


3 

v 


where 


2  4  6 

V  V  V  V 

IT  60  "5040  +  272160 


(A.  3) 


C2n  _  1 . 4380n 

T  T  (A.  4) 


This  series  converges  for  v<  2»rbut  convergence  is  slow 
when  v  >  1.  Defining 
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-  .  V2  +  V1  .  .  V2  -V1 
V  =1  - - -  ;  A  v  =  - - - 


(A.  5) 


Equation  A.  3  was  transformed  to 


V  = 


15 


ir 


( 


Aviv2 


-3  -  4 

+  V 


-  6 

_  _v _ 

12  "  720 


+  (A  v)’ 


+  (A  v)' 


-  2 
v _ 

24 


(-J - 

\  12  8 

/  _i__  _ v.2  +  \ 

Y  960  3840  *•**  j 


-4 

v 

576 


^  •  •  •  « 


+  0  [  (Av) 


(A.  6) 


and  used  to  calculate  t)  when  v  £  1.  Inclusion  of  the  terms 
shown  through  order  (A  v)^  results  in  a  value  of  Tf  accurate  to 
±0.1%. 


2.  While  D  itself  cannot  be  expressed  in  closed  form,  its 
derivatives  can- be.  Since  ?7  is  the  difference  between  two 
values  of  D,  the  zero  order  term  of  a  Taylor  expansion  of  TJ 
vanishes.  Expansion  about  v  eliminates  all  even  order 
derivatives,  thus 


r?  =  [  (v2-v)  -  (v2  -v}  ] 


a  D 
a  v 


^v  =  v 


3  ! 


(Vj  -V)3  -  (v2-v)3 


a3D 


a  v~ 


v=v 


(A.  7) 


Since 

a  p 
a  v 


15 


(e  -1) 


(A.  8) 
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the  energy  fraction  becomes 


+  {A  v* 


hibvll 


(6+v^+6v~I2)  eV+  (v^-6v-f6)  e^V 


Retaining  only  the  term  of  order  A  v,  Eq.  (A.  9)  was  used 
to  calculate  T}  to  ±  0.  6%  when  v  >  1.  The  curve  v  =  1  on  Fig. 
3.  1  separates  the  region  of  application  of  the  two  methods. 
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APPENDIX  B 
PROGRAM  ARCRAD 


B.l  PURPOSE 

The  purpose  of  this  program  is  the  calculation  of  radiation 
heat  flux  from  a  cylindrical  air  arc  based  on  spectral  emiBsion- 
absorption  characteristics  of  air  and  a  measured  or  assumed  tem¬ 
perature  profile. 
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B.2  GLOSSARY  OF  TERMS 

Fortran  Name  Formula  Name  Description 


ABE  MIS  (2  2) 

AINC 

ALPHA(IO)  x 

AMAX 

AMIN 

AKDR(45) 

ANGADiD  x/  2 

ANGLE  (20) 

B 

BCDREC(14) 

BUFF2(525) 

BUFFT(86) 

C 


Intermediate  array  used  in 
calculation  of  EOUT. 

Decrease  in  wave  number. 

Angles  used  for  absorption 
lengths  for  calculation 
EXTRMN  and  EXTRMX  and 
EIN. 

Running  sum  of  maximum 
radiation  heat  flux  from  all 
annuli  of  cylinder  (excluding 
core)  -  final  sum  corresponds 
to  only  one  wave  number  WN(I). 

Running  sum  of  minimum 
radiation  heat  flux. 

Components  of  linear 
absorption  path. 

Angle  difference  in  radians  of 
elements  of  ANGLE  array. 

Angles  used  for  absorption 
lengths  and  transmissivity  for 
interior  radiation. 

Intermediate  value  in  calculation 
of  RI  array. 

Storage  locations  for  reading 
BCD  input. 

Buffer  for  labeled  and  blocked 
tape  on  A6. 

1-0  buffer  storage  total. 

Intermediate  value  in  calculation 
RI  array  RI(NQ,  1 ). 
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Fortran  Name 

Formula  Name 

Description 

ENTRMN 

Ti,  LL 

Minimum  transmissivity  for 
interior  radiation. 

ENTRMX 

Ti,  UL 

Maximum  transmissivity  for 
interior  radiation. 

EOUT 

€e 

Emissivity  exterior. 

EFSMN 

Minimum  interior  emissivity. 

EPS  MX 

Maximum  interior  emissivity. 

ETA(23,  121) 

V 

Fraction  of  total  energy  radiated 
by  black  body  at  annulus  temper¬ 
ature  and  wave  number. 

ETASM(23) 

• 

ETA  of  each  temperature  summed 
over  all  wave  numbers. 

EXRQSN(22,  10) 

Intermediate  array  used  in 
calculation  of  R  array. 

EXSNRQ(22,  10) 

Intermediate  array  used  in 
calculation  of  R  array. 

EXTRMN 

Te,  LL 

Minimum  transmissivity  of 
exterior  radiation. 

EXTRMX 

*e,  UL 

Maximum  transmissivity  of 
exterior  radiation. 

FKR{20) 

Empirical  approximation  of  one 
of  the  integrals  used  in  calculation 
of  transmissivity  and  interior 
emissivity. 

FMT(12) 

Storage  for  input  Hollerith 
identification. 

I 

Index  for  reading  in  data  on 
tape  A2,  for  printing  output  on 
tapes  A3  and  A6,  and  indexing 
wave  numbers. 
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Fortran  Name 

11 

12 

IAL 

IOPT  1 

IP  1 
IPST 

J 

JAL 

K 

LAL 

MAX- 

N 

NA 

NAH 


Index  used  in  reordering  the 
part  of  TRYKDR  array  - 
associated  with  2  consecutive 
ALPHA  -  in  ascending  order. 


Index  used  in  reordering  part 
of  TRYKDR  array  in  ascending 
order. 


Subscript  used  in  ordering 
linear  absorption  path  for 
interior  radiation. 


Option  for  printing  on  A6  - 
save  tape  used  in  printer. 


Lower  limit  of  12  index. 


Adjustment  value  used  in 
selecting  linear  absorption 
coefficients  from  proper  input 
records. 

Subscript  used  to  find  annulus 
position  -  used  in  determining 
AKDR  array  for  interior 
radiation. 

Upper  limit  for  II  limit. 

Index  used  for  reading  record 
of  files  of  input  tape. 

Index  for  reordering  entire 
TRYKDR  array. 

Lower  limit  of  II  index. 

Index  for  annuli. 

Index  for  ALPHA  array. 
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Fortran  Name  Formula  Name 

NAL 

N ALPHA 

NANGLE 

NAVE 

NCASE 

ND1V 

NFAC 

NFILE 

NH 

NJ 

NPP 

NQ 

NRAVE  s 


Description 

Upper  limit  of  12  index. 

Number  of  members  in 
ALPHA  array. 

Number  of  members  in 
ANGLE  array. 

Index  for  ANGLE  array. 

Number  of  cases  of  data  to  be 
used  for  a  computer  run. 

Number  of  sectors  of  half 
cylinder. 

NTOUCH+1  -  the  first  annulus 
counting  from  core  outward 
that  the  linear  absorption 
path  of  interior  radiation  inter¬ 
sects. 

Counter  for  file  number  used 
during  read  in  process. 

Value  used  as  part  of  subscript 
in  calculation  of  EXTRMN. 

Subscript  used  in  ordering  and 
summing  the  components  of 
linear  absorption  path  of 
interior  radiation. 

Upper  limit  of  K  index. 

Index  used  in  calculation  of  R 
along  linear  absorption  path  of 
both  interior  and  exterior 
radiation  -  also  used  in  interior 
radiation. 

Number  of  divisions  of  the 
sectors  of  half  cylinder. 
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Fortran  Name 

Formula  Name 

Description 

C2T(23) 

tf/T 

Intermediate  array  used  in 
calculation  of  V, 

CORMN  (121) 

Minimum  radiation  heat  flux  of 
core  -  corresponds  to  one  wave 
numbe  r. 

COR  MX  (121) 

Maximum  radiation  heat  flux 
of  core. 

D 

Intermediate  value  used  in 
calculation  of  RI  array-RI(NQ,  2) 

DALPHA 

n{  2s 

Angle  difference  in  degrees  of  the 
elements  of  ALPHA  array. 

DELTA 

s 

Radius  of  arc  column. 

DK 

Absorption  length  along  diameter 
of  arc  column. 

DR(23) 

Annulus  thickness  -  outermost 
inward. 

DRK 

Absorption  length  at  0°  used 
in  exterior  radiation. 

DSINA(20) 

Sine  of  angles  of  ALPHA  array. 

DV(23) 

Intermediate  value  in  calculation 
of  ETA  array. 

EIN 

Emissivity  inward. 

EKDR(IO) 

Absorption  length  used  in 
calculation  of  EIN. 

ENSNA(20) 

Sine  of  angles  of  ANGLE  array. 

ENSNRQ(23,  20) 

Intermediate  array  used  in 
calculation  of  RI  array. 

ENRQSN(23,  20) 

Intermediate  array  used  in 
calculation  of  RI  array. 
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Fortran  Name  Formula  Name 

NR 

NTEST 

NTOUCH 

NT  m 

NTR 

NUMDR 

NWB1N 

NWBCD 

NW 

OPT 

OPTCOR(121) 

QMAX{121)  W 

UL 


De  scription 

NT+1  number  of  annuli, 
plus  one. 

Number  of  annuli  {depending 
on  position  of  radiating  annulus) 
including  core  and  going  out¬ 
ward  which  may  be  intersected 
in  two  positions  -  used  in 
calculation  of  RI  for  interior 
linear  absorption  path. 

Last  annulus  (considering  core 
as  annuli  one)  not  intersected 
by  interior  linear  absorption 
path. 

Number  of  annuli  of  constant 
temperature  of  the  cylinder, 

NT-1  number  of  annuli  minus 
one. 

Number  of  components  of  interior 
linear  absorption  path. 

Control  parameter  =  0  to  read 
in  entire  binary  word. 

Control  parameter  =  0  to  read 
in  entire  BCD  record. 

Number  of  wave  numbers. 

Running  sum  of  radiation  heat 
flux  from  all  annuli  (not  including 
core)  -  assuming  gas  non-self 
absorbing. 

Radiation  from  core  assuming 
gas  non- self  absorbing. 

Radiation  heat  flux  per  wave 
number  at  boundary  of  cylinder 
considering  maximum  trans¬ 
missivity. 
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Fortran  Name 


Formula  Name  Description 


QMIN{121) 

QOPT(121) 

QOPTSM 

QSUMMN 

QSUMMX 

RAD(23) 


W 


LL 


W. 


OP 


r 


Radiation  heat  fLux  per  wave 
number  at  boundary  of  cylinder 

-  considering  minimum  trans- 
mis  sivity. 

Radiation  heat  flux  per  wave 
number  at  boundary  of  cylinder 

-  assuming  gas  non-self 
absorbing. 

Radiant  heat  flux  QOPT(I) 
summed  over  all  wave  numbers. 

Minimum  radiant  heat  flux 
QMIN{I)  summed  over  all  wave 
numbers. 

Maximum  radiant  heat  flux 
QMAX(l)  summed  over  all  wave 
numbers. 

Array  of  radii  ordered  from 
outermost  annulus  in  to  core. 


RADQ(23) 


Array  of  radii  order  from  core 
to  outermost  annulus. 


RAVE 

RADLPH 

RDE  LTA(23) 

REC(700) 
RHO  (23) 
RI(23,  2) 


s 


log  P/Po 

R+ 

R' 


Floating  NRAVE. 

Increment  in  radians  of  elements 
of  ALPHA  array. 

Array  of  radii  of  annuli  over 
radius  of  cylinder  -  from  outer¬ 
most  annulus  in  to  core. 

Storage  for  reading  in 
binary  words. 

Log  of  density  of  annuli  from 
outermost  annuli  in  to  core. 

Sector  radii  array  along  path 
of  interior  radiation. 
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Fortran  Name 
R(22) 

SIGMA 

SIN  A 
SKDR 

SMCRMN 

SMCRMX 

SMOPTC 

SUM 

TEMP  (23) 

TEMPL 

TEMPRW{23) 

TEMPS 

V 


Description 

Sector  radii  along  path  of  exterior 
radiation. 

Intermediate  value  in  calculation 
of  RI  array. 

Sine  of  ALPHA  array. 

KR  Dimensionals  absorption  length 

used  in  calculation  of  trans- 
mis  sivity. 

Sum  over  all  wave  numbers 
of  minimum  radiant  heat  flux 
core  (CORMN(I)). 

Sum  over  all  wave  numbers  of 
maximum  radiant  heat  flux 
from  core  (CORMX(I)). 

Sum  over  all  wave  numbers  of 
radiant  heat  flux  from  core 
(OPTCOR(I))  -  assuming  gas 
is  non  self-absorbing. 

Intermediate  value  used  in 
summing, 

T  Temperature  array  -  corresponding 

to  outermost  annulus  in  to  core. 

Temporary  storage  used  in  re¬ 
ordering  of  TRYKDR  array. 

Intermediate  array  used  in 
calculation  -  o  T^. 

Temporary  storage  used  in 
reordering  TRYKDR  array. 

v  Parameter  using  Plank  radiation 

function  for  the  calculation 
of  E  TA  array. 


Formula  Name 
R 
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Fortran  Name 

Formula  Name 

Description 

V2 

Second  power  of  V. 

V3 

Third  power  of  V. 

V4 

Fourth  power  of  V. 

V6 

Sixth  power  of  V. 

WBB 

Intermediate  value  in 
calculation  of  WIN  and  WOUT. 

WIN 

Radiation  inward  -  gas  non 
self-ab  sorbing. 

WINMN 

Radiation  inward  -  considering 
minimum  transmissivity. 

WINMX 

Radiation  inward  -  considering 
maximum  transmissivity. 

WNA 

Highest  wave  number  plus  AINC. 

WN(121) 

X 

Wave  number  array. 

WOUT 

Radiation  outward  -  gas  non 
self-absorbing. 

WOUTMN 

Radiation  outward  considering 
minimum  transmissivity. 

WOUTMX 

Radiation  outward  considering 
maximum  transmissivity. 

X 

Intermediate  value  in  calculation 
of  ETA. 

Y(4) 

Array  used  in  printing  on  save 
tape  A6. 

Z 

Intermediate  value  in  calculation 
of  ETA. 
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B.3  INITIALIZATION  AND  INPUT,  PROGRAM  NOTES 
B.3.1  Initialization 

1.  General  Initialization 

The  Initial  sequence  sets  aside  86  locations  in  the  BlJFFT 
array  for  the  1-0  buffer.  Input  requires  39  locations- 
output  47  locations. 

2.  Input  Initialization 

Call  ABUF1  sets  up  input  buffer  on  A5.  To  read  BCD, 
decimal  tape  number  is  645.  To  read  BIN,  decimal  tape 
number  is  661.  Logical  tape  number  is  15. 

3.  Output  Initialization 

BUFFS  sets  up  output  buffer  on  A6  to  write  labeled  and 
blocked  tape.  To  read  BIN,  decimal  tape  number  is  662. 
Logical  tape  number  is  16. 

4.  Header  for  File  of  LBT  on  A6 

Call  LBID,  specifies  a  BCD  file  name  for  label  of  a 
labeled  and  blocker  tape  to  be  written  on  A6.  Call  LBID 
for  each  file. 

5.  Identification  of  Case 

Call  SPGHDR  (FMT):  page  header  is  written  with  12  BCD 
words  of  array  FMT.  Call  SPGHDR  for  every  case. 

6.  Error  Initialization 

ERTRP  allows  printout  and  transfer  to  next  case 

B.3.2  Input 

1.  Tape 

Binary  tape  data  is  input  for  temperature  3000°K-25000°K. 
This  is  a  tape  with  23  files,  one  file  per  temperature.  A 
description  of  tape  will  be  given  in  B.  3.  3. 
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2.  Cards 

a.  number  of  cases  (fixed  point) 

b.  Cards  per  case 

1)  identification  card 

2)  Radius  of  cylinder,  decrease  in  wave  number, 
highest  wave  number  -  AINC,  angle  difference 

in  degrees  of  array  ALPHA,  divisions  of  original 
sectors  (ALPHA  array)  of  half  cylinder  (floating 
point) 

3)  number  of  wave  numbers,  number  of  annuli  or 
temperatures,  option  for  writing  output  tape,  NT-1* 
NT+1,  number  of  angles  in  ALPHA  array,  divisions 
of  original  sectors  of  half  cylinder  (fixed  point) 

4)  temperature  array:  temperatures  of  annuli  from 
outermost  temperatures  must  be  multiples  of 
thousands  in  ascending  order  into  core  (floating 
point). 

5)  radius  array:  radii  of  annuli  from  outermost  into 
core+1  (floating  point) 

6)  radius  array:  radii  of  annuli  from  core  to  outer¬ 
most  (floating  point) 

7)  log  of  density  ratio  array  -  log  of  density  ratio  of 
annuli  from  outermost  into  core  (floating  point) 

B.3.3  Read  in  Process 

Read  in  linear  absorption  coefficients  for  each  case*  inter¬ 
polating  as  necessary. 

1.  General 

A  table  of  absorption  coefficient  data  is  interpolated  with 
density  selected  from  input  tape  and  set  up  as  a  function 
of  temperature  and  wave  number,  f  (TEMP*  WN) 
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2.  Format 

Data  is  stored  in  files.  Each  file  corresponds  to  a  tem¬ 
perature.  File  temperatures  range  from  3000°K  to 
25000°K  in  1000°K  increments. 

The  first  file  (3000°K)  is  made  up  of  5  records.  Record 
1  is  a  BCD  label.  Record  2  is  main  table  heading.  The 
main  table  headings  are:  NT,  TEMP(l),  TEMP(2)„.. 
TEMP{23).  Record  3  contains  the  number  of  words  in 
the  data  subtable,  plus  heading:  NW,  ND,  FMIC(l) 

FMIC(2). .  .  FMIC(NW)  DENS{1),  DENS(2). . .  DENS(ND). 

Note:  NT  is  number  of  files 

NW  is  number  of  wave  numbers  (21) 

ND  is  the  number  of  densities  (5) 

FMIC(l)  is  micron  value  of  WN(1)  =  10/WN(1) 

FMIC(NW)  is  micron  value  of  WN(NW)  =  10/ WN(NW) 

DENS(l)  =  -3  (RHO/RHOO) 

DENS(ND)  =  1  (RHO/RHOO) 

Record  4  consists  of  one  word,  the  .temperature  of  the  file. 
The  last  record,  record  5,  is  the  data  subtable  and  con¬ 
tains  linear  absorption  coefficients  for  the  temperature  of 
the  file  for  all  wave  numbers  for  each  of  the  5  log  density- 
ratios.  The  remaining  22  files  are  composed  of  3  records, 
corresponding  to  the  format  of  the  last  3  records  of  file  1: 
a  subtable  record,  a  temperature  record,  and  a  data  re¬ 
cord. 

Temperatures  (multiples  of  1000°K)  in  array  TEMP  are 
stored  in  ascending  order  as  are  the  file  temperatures. 
Reading  of  files  is  controlled  by  DO  N=l,  NT  (number  of 
temperatures).  Records  of  files  are  read  up  to  and  in¬ 
cluding  temperature  record.  If  the  first  and  only  word 
of  the  temperature  record  is  not  equal  to  the  indexed 
TEMP(N),  the  file  is  closed  and  the  next  file  is  then  read 
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in  a  similar  manner.  If  TEMP(N)  and  the  one  word  of  the 
temperature  record  are  equal,  the  following  record  is 
read  and  processed  for  all  wave  numbers  (61500-1500) 
to  obtain  the  absorption  coefficients.  PROP(N,  l)  through 
PROP(N,  121)  for  TEMP(N)  with  its  associated  log  of 
density  ratio  (RHO(N).  The  file  is  closed,  the  index  for 
TEMP(N)  is  increased  by  1  and  the  reading  process 
continued  to  obrain  linear  absorption  coefficients  for  all 
the  indexed  temperatures.  The  first  file  is  read  as  a 
special  case.  The  BCD  record  is  read  by  a  special  in¬ 
struction  and  then  the  3  following  records  must  be  read 
to  include  the  temperature  record.  The  remaining  22 
files  are  similar  in  read  instructions;  2  records  are  read, 
the  second  being  temperature  record. 


Records  are  all  read  in  REC(700)  which  is  large  enough  to 
include  the  largest  record  (the  last)  of  the  files. 


The  data  for  TEMP(N)  is  stored  as  a  function  of  density 
and  wavenumber  in  the  following  order:  f(d.,NW.) 
f(d2,NW1)....f(d5,NW1),  f{d1(NW2),  f(d2>lW2)/... 

fldj.NW^),  . .  .f(dg,  NW121).  Where  d  is  density  and  NW 
is  wave  number.  Thus  the  data  record  (last  record  of 
each  file)  can  be  treated  as  121  groups  of  5.  Each  group 
corresponding  to  one  wave  number.  The  5  linear  absorption 
coefficients  in  each  group  are  arranged  in  ascending 
order  of  log  density  ratio  -3,  -2,  -1,  0,  1. 


If  the  RHO(N)  corresponding  to  TEMP(N)  is  an  integer 
value  between  -3  and  1,  the  value  IPST,  which  is  the  sub¬ 
script  of  array  REC  is  assigned  as  follows: 


IPST  =  1  if  RHO(N)  =  -3 
IPST  =  2  if  RHO(N)  =  -2 
IPST  =  3  if  RHO(N)  =  -1 
IPST  =  4  if  RHO(N)  =  0 
IPST  =  5  if  RHO(N)  =  -1 
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Control  is  then  transferred  to  DO  loop  -  DO  1  =  1, 

NW(121).  The  value  REC(IPST)  is  read  directly  for 
linear  absorption  coefficient  for  first  wave  number. 

IPST  is  increased  by  5  in  last  step  of  DO  loop.  When 
the  DO  is  satisfied  the  linear  absorption  coefficients 
PROP(N,  1)  through  PROP(N,  121)  are  obtained  for 
TEMP(N)  with  its  corresponding  RHO(N). 

If  RHO(N)  is  not  an  integer  value  but  within  the 
integer  range,  the  linear  absorption  coefficient  is  ob¬ 
tained  by  logarithmic  interpolation  between  tabulated 
values.  IPST  is  set  corresponding  to  the  next  lower 
RHO(N),  e.g.  RHO(N)  =  -2.5,  1PST  =  1>PCNT  is  the 
fractional  difference  between  tabulated  values,  REC(IPST) 
and  REC  (IPST  4-1).  PC  NT  is  equal  to  RHO(N)  minus 
integer  values  less  than  RHO(N).  Control  is  then  trans¬ 
ferred  to  a  DO  loop  -  DO  1=1,  NW(121).  Linear  ab¬ 
sorption  coefficient  for  first  wave  numbers  equals  antilog 
of  PC  NT  [log  [REC  (IPST  +  1)]  -  log  [REC  (IPST)]  }  + 
log  l  REC  (IPST)  ].  IPST  i9  increased  by  5  in  last  step 
of  DO  loop.  When  the  DO  is  satisfied  the  linear  absorp¬ 
tion  coefficients  to  PROP(N,  1)  through  PROP(N,  121) 
are  obtained  for  TEMP(N)  with  its  corresponding  RHO(N). 

If  RHO(N)  is  greater  than  1,  PCNT  is  calculated  as  pre¬ 
viously  stated,  and  IPST  is  set  equal  to  5.  Absorption 
coefficient  is  logarithmically  extrapolated  based  on  values 
of  REC  (IPST)  and  REC  (IPST  -  1). 

The  log  difference  of  REC(IPST)  and  REC(IPST-l)  is 
added  to  REC(IPST)  and  an  antilog  is  taken  to  obtain  in 
DO  loop  similiar  to  the  one  previously  described  PROP 
(N,  l)  through  PROP(N,  121)  for  TEMP(N)  with  its  corre¬ 
sponding  RHO(N). 

This  concludes  the  read  in  part  of  the  program.  Notes  related 
to  the  remainder  of  the  program  are  incorporated  in  the  block  diagram, 
section  B.  4. 
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1 - 

calculate  intermediate  arrays  used  in  calculation  of  absorption 

NAVE  =2, 

lengths  of  exterior  radiation 

NANGLE 

ENRQSN(NQ.NAVE) 

1 

ENSNRQ(NQ,  NAVE ) 

- 
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96 


calculate  V  =  WNJI)*C2T{N)  Calculations  are  limited  to 
the  wavenumber  range  which  includes  all  but  0. 1%  of 
black  body  energy  at  each  frequency  extreme  by  trans¬ 
ferring  to  end  of  do  loop  when  v  falls  outside  the  range 
(.278, 13).  Whether  V*  1  or  >  l  determines  the  equation 
for  calculating  ETA(N,  I),  the  fraction  of  total  energy  flux 
radiated  by  black  body  at  a  given  temperature  within  a 
wave  number  interval. 


}C ALL  ERROI 


N<JNT 


256  ERASE  SK DR 


CALCULATION  OF  ABSORPTION 
LENGTH  FOR  EXTERIOR  TRANS¬ 
MISSIVITY  AND  INTERIOR 
EMISSIVITY 


calculate  sector  radii  R(NQ)= 


CALL  ERROR 


AEDC-TR-65-1  1 


N> 


CALL  ERROR}? - WQri - 


NQ-1 


NQ=1,N 


NQ>  1 


NAt2, 

NALPHA 


ich| 

d 

to 

us 


260  (  AKDR(NQ-i)=(R<NQ-l)-R(NQ)?-yROP(NQ-l.I)| 


CALL  ERRORS — NQ} 


NQCN- 


NQ 


N 


NA  designates  angle 
of  wl| 
is 

ward 
annu 
radii  tied  from  the 
Nth  i  nnulua. 


subscript  NQ 
<  minting  in- 
lenote  the 
Lb  sorbing 


AKDR(N)=R(N)=fPROP(N,I) 

EKDR{NA)=AKDR{"NJ 

EKDR  are  absorption  lengths  associated  with 
TEMP(N)  for  aU  ALPHA |N A)  for  the  calculation 
of  emissivity  inward 

AKDR  are  components  of  the  absorption  path 


-2W - .| CONTINUE  jr 


NQ*1P  N 

L_ 


SKDR(NA)=SKDR(NA)+AKDR(NQ) 

SKDR  are  absorption  length  associated  with 
TEMP(N+1)  for  all  ALFHA(NA)  for  calculation 
of  transmissivity  (exterior} 


SKDR<1)=DRK 
absorption  length 
at  ALPHA  (1) 
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CALCULATION  OF 
INTERIOR  EMISSIVITY 


FKR(NA)  is  an  empirical 
approximation  to  the  trans¬ 
missivity  integral;  the  aolutijo: 
of  FKR{NA)  varies  if  EKDR 
1 1.  or  >  16. 


(NA) 


Interior  emissivity  of  a  single 
annulus  la  identical  to  transmissivity 
of  that  annulus  alone  relative  to  the 
next  interior  annulus,  r 


ERASE  EPSMN,  EPS  MX 


NA=2, 
NA,LPHA 


-4- 


310 


Transmissivity  is  sum  of 
products  of  2  integrals,  one 
being  DSINA(NA)  the  other 
approximated  by  FKR(NA) 

The  lower  limit  of  trans¬ 
missivity  corresponds  to 
sector-segment  coincidence  at 
NA=2,  ■  •  •  NALPHA  -  upper 
limit  at  NA=  1 .  •  •  •  NALPHA  - 1 


Prior  to  the  Do  statement  EXTRMN  and 
EXTRMX  are  Bet  equal  to  1,  and  thiB  is  the 
value  used  when  the  Do  loop  index  has  its 
initial  value  (N=l),  Later  in  the  loop 
EXTRMN  and  EXTRMX  are  calculated 


EPSMN=EPSMN+FKR(NA)*LSINA(NA) 

E  PS  MX=E  PS  MX+FKR  (N  A- 1 ) <DS  IN  A(  NA ) 


emissivity  exterior  to  annulus 
EOUT=l.  -EXPF(-”*) 


interior  emissivity  is  the  average 
of  upper  and  lower  limit 
EIN=1, -(EPSMN+EPSMX>/2r 

- T - 

_ if— _ _ _ 

WBB  =E  TA(N,  I)*TEMPRW(N) 
WIN=WBB<ErN*iRDELTA(N) 
WOUT=WBB<£OUT*RDELTA(N} 
WOUTMN=WOUT«XTRMN 
WOUTMX=WOUT«EXTRMX 


99 


AE  DC-TR-65-1 1 


CALCULATION  OF 
ABSORPTION  LENGTH 
FOR  INTERIOR  TRANS¬ 
MISSIVITY 


■*1 

ERASE  TRY 
array  of  abs 
N)  associate 

: _ 

KDR 

orption  length  (for  annulus 
d  with  ANGLE  array 

NTEST=NT-N-1 

NTEST=#  of  sector  radii  lees  than 


RAD(N  -11 ) 

annulus  NTEST+1  and  all  succeeding 
annuli  can  be  intersected  by 
radiation  path  only  once 


Because  absorption  coefficients  are  not 
monotonically  increasing  function  of  temperature, 
the  extremals  of  absorption  length  within  a  sector 
do  not  always  correspond  to  the  value  at  sector 
boundaries  ALPHA(NA)  and  ALPHA(NAtl), 

An  absorption  coefficient  is  also  calculated  within 
each  sector  and  the  smallest  and  largest  of  the  3 
hre  used  in  computing  maximum  and  minimum 
transmissivity  respectively.  The  ANGLE  array 
includes  the  angles  of  the  ALPHA  array  and 
those  angles  midway  between  successive 
elements  of  the  ALPHA  array. 
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NAL=NAL-rNRAVE 

NRAVE  read  in  as  2  (division  of  original  angle 
difference  used  for  exterior  radiation) 
TEMPL=TRYKDR(NAL) 

store  the  absorption  path  corresponding  to 
ALPHA  (2)  to  ALPHA(NALPHA)-upper  boundary 
of  sector 
JAL=NAL-1 

upper  limit  of  index  reordering 


IP  1=11  -1 


DO  1 1  =MAL,  JAL 


lYKDR(Il)  TRYKDRCE] 


DO  IZ=IP1,  NAL 


When  the  DO  has  been 
satisfied,  TRYKDR(MAL) 
to  TRYKDR(NAL)  has  been 
rearranged  so  that  the 
smallest  absorption  path 
is  in  TRYKDR(MAL). 


TEMPS=TRYKDR(I2) 
TRYKDRd  I )  =TRYKDR(E) 
TR  YKDR  (12 ) = TE  MPS 


(first  element  is  compared  with 
second  and  third  and  inter¬ 
changed  when  necessary) 


CONTINUE 


This  DO  governs  the  reordering 
from  TRYKDR(MAL)  to 
TRYKDR(NAL)  in  ascending 
order:  TR  YKDR  (NAL)  is 
largest  absorption  path. 


First  halt  of  SKDR  array  SKDR(l)  to  SKDR(KDIV)  will 
be  stored  with  smallest  of  interior  absorption  path 
between  ALPHA(NA)  and  ALPKA(NA+1);  the  remaining 
half  SKDR(NDrV+l)  to  SKDR(NDIV+NDlV)with 
largest  interior  absotption  path  between  ALPHA(NA) 
and  ALPHA(NA+1) 
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'  I - 

NA=2, 

NALPHA 


inte  riorJtransmiasivity  to  be  used  with  TEMP(N) 
NAH=NA+NH 

NAH  is  subscript  used  to  pick  the 
largest  interior  absorption  coefficient 


of  each  sector 

ENTRMN=ENTRMN+FKR(NAH)*DSINA(NA) 

ENTRMX=ENTRMX+FKR|NA-l)’iDSINA(NA) 


radiant  heat  flux  WINMN  (minimum) 
and  WINMX  (maximum)  for  annulus 
N  considers  only  transmissivity 
and  emlssivity  interior 

radiant  flux  density  interior 
and  exterior  to  annulus  is  kept 
as  a  running  sum  over  all  the 
annuli  excluding  core -minimum 
and  maximum  transmissivity 
considered  separately  -  case 
of  no  absorption  is  considered 
In  a  similar  manner 


WENMN^WIN^ENTRMN 

WINMX=WIN*£NTRMX 


I 

_ _ _  - 

AMIN = AMIN+  WIN  MN+WOUTMN 

AMAX=AMAX+WINMX+WOUTMX 

OPT=OPT+WIN+WOUT 


650 


CONTINUE 


stored  in 
array  as 
function  of 
wave  number 


700 


option  for  writing 
data  for  case  on  tape  A6 
saved  and  used  for 
plotting 


OPTCOR(I)=-*  • 

radiation  from  core-no  absorption 
CORMN(I)=OPTCOR(l)«XTRMN 

radiation  from  core-min,  trans¬ 
missivity 

CORMX(I)=OPTCOR(I)«£XTRMX 

radiation  from  core  -  max.  trans- 
mis  sivity 

QMIN(I)=AMIN+CORMN(I) 

minimum  radiation  reaching  boundary 
QMAX([)=AMAX+CORMX(I) 

maximum  radiation  reaching  boundary 


QOPT(1)=OPT+OFTCOR(I) 

radiation  reaching  boundary  not 
considering  absorptivity 


i 

i 
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YES 


arrange  data  in  special 
format  to  be  written  on  tape 

wave  number  and  radiation 
heat  flux  as  function  of  wave 
number  written  in  separate 
records 


_ it _ 

ERASE  Y 

Yjl)=WN(I) 

Y(2)=QMAX(I) 

Y(3)=QMIN(I) 

Y(4)=Q0PT(I| 

CALL  RITE  (662,  Y,  4) 


go  to  EOF  on  output  tape 
after  writing  data  for  each 
case 


± _ _ 

CALL  CLOSE(662)  j 


WRITE  VALUES  IN  BCD  ON  A3 


write 


(WOT  10,  1  \ 
WOT  10,  B  'l 
WOT  10.  7  / 


wI  nT  JffOT  10,  9,  TEMP(I),  RHO(I),RAD(I)j 


summation  of 
radiation  minimum, 
maximum  and 
optically  thin  - 
summation  of 
radiation  from  core 
in  same  manner 


write  /WOT  lO.TN 
headers  \WQT  10,4  J 


1 


=1 ,( NW  [WO T  1 0,  5,  WN(I),QMAX{I),QMIN(l).QOPT(I); 


ERASE  QSUMMN.QSUMMX.QOPTSM, 
SMCRMN.SMCRMX.SMOFTC 


QSUMMN=QSUMMN+QMIN(I) 

QSUMMX=QSUMMX+OMAX(I> 

Q  OPTSM=QOPT5M+CX)PTaj 

SMCRMN=SMCRMN+CORMN(I) 

SMCRMX=SMCRMX+CORMX(I) 

;  SMOPTC=SMOPTC+OPTCOR(I) 

- | - 

I 

| 

i 


i 

i 
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fraction  of  total 
energy  flux 
radiation  by 
black  body 
within  wave 
numbe  r 
interval  WN(I) 
to  WN(NW) 


(  WOT  10,  13,  SMCRMX, 
\  WOT  10,  6,  QSUMMX.C 


SMCRMN.SMOPTC 
QSUMMN,  QOPTSM 


-4 

r 

ERASE  ETASM 

_ _ _ - _ b 

t _  ... 

, - 1 - 

=1,  NT  1=1,  NW 

ETASM(N)=ETA(N,I)+ETASM(N)  | 

"4 

L _ 

) 


WOT  10, 10  ) 


-I 


CONTINUE 


H 

MZE  661,,  32  384 

K 

MZE  662,  ,  32384 

rewind  and  unload  tape 

CALL  \ 
EXIT  ■ 
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B.5  PROGRAM  LISTING 

STUDY  nF  RADIATION  HEAT  FLUX  FROM  HIGH  PRESSURE  AIR  AKCS  AT  LCM  F09/ L 1/6A 


c 

c 

COMMON  A0EKIS(22),AINC,AKDR( A5I, ALPHAI  10  I (ANGLE <201 
COMMON  DCOHECl 25) 

COMMCN  fiUFF2(525)»8UFFT(86),  CORHN  (121)»CORMX(121)  ,C2T(23) 

COMMON  DALPHA, DELTA,  CPt- l 

COMMCN  DR  I  231 ,DS [NAI 10 > , C V < 2 3 J , EKDR l 10 1 , ENSNA 1 2 0 ) 

COMMCN  CNSNRQt  23, 20  , ENRCSN(  23,  20  I,  ETA  123, 1 21 T 
COMMON  E  T  ASM ( 23  I 

COMMON  EXSNRQI22,20),EXRCSNI  22 , 2C ) , FKR < 20 ) ,  FM  T  (  12  >  ,  I T  {  2  3 ) 

COMMON  I  CPF  1 ,  NALPFA , NO [ V 

COMMON  iMRAVE.NTR.NT.NH.OPTCURl  121) 

COMMON  PRCP( 2  3, 1211 ,CM«X( 121 1 , CM  INI  121) 

COMMON  OCPI ( 121 ) »RAOC(23)»ftADl23l«RAVE,ROALPH,ROELTA(23) ,REC(700) 
COMMON  RHOI 231 

COMMON  »II 23,2) ,R(22).SINAI20),SKDRI20) .TEMP (23) , TEMPRW (23) 

COMMON  TRYKDRl 20 ) , WN A , WN ( 1 2 1 1 , V <  A ) 

1  FORMAT ( 1 H 1 ) 

2  FORMAT ( 12A6 ) 

A  FURMAT  IIHC3X0HWAVE  NO. , 1 LX4HCMAX , 13XAH0M1N, 13X4H0GPTI 
5  FORMA T ( IH  M 1PE1A.7.3X)  1 

t  FURMAT( 1HQ3X0HSQMAX  =  1PE 14. 7, 3X 8HSSM1N  =  1PE 1 A . 7 , 3 X8HSQCPT  -  1PEI 
1 A  .  7  ) 

7  FURMAT  [  LHC2X11HTEMPETATURE,  5X12HLOG  R-HU/RHCO,  8X6HRADIU5) 

8  FORMAT!  LH03 0 X 19HT EMPERA TLRE  PROFILE) 

S  FCRMAT ( IH  3 (  IPEIA.7, 3X 1  ) 

1C  FORMAT  ( 1HC2X1 IHTEUPERA TORE, 8X6HETASLM) 

11  FORMAT  l  IH  2UPE1A.7.3X)  ) 

13  FORMAT  < IHCIXIOHSQCORMX  =  IP E 14 . 7 , l X 10HSQCCRMN  =  1PE1A.7, 

11X1 1FSCCCRC.PT  =  1PE1A.71 
C  INITIALISE  FOR  INPUT  OUTPUT 

5C  IF(  INI T(6A2,6A3,BUFFT (8  6) ,  39, A7 1)60,55,60 
55  CALL  DUMP 

C  INITIALIZE  FOR  ERROR  ROUTINE 
60  IF  IERTKPI5) 170,726,726 
C  IM TIALISE  FOR  BUFFERS 

C  AEUFl  IS  FOR  INPUT  ON  A5 
TC  CALL  A8UFI 

RJT  2,2.  NCASE 

IF  f8UFFS(662,0UFF2( 525 1,525,0, 1, 1, l )) Tl, 55,71 
C  READ  IDENTIFICATION  CARD 

C  REAC  IN  CARDS  CN  A2 

71  RIP12,2,(FMTt J), 1=1, 12) 

CALL  SPGHDR  I FMT ) 

C  INPUT 

R IT  2, 2, CEL T A, A  INC, WN A, DALPHA, RAVE 

RIT  ^,2,  Nh.NT, IOP T 1 ,NT R , NR, NALPHA , NR AVE , ND I  V 

RIT  2,2,  (TEMPI  I), 1  =  1, NT) 

RIT  2,2,  I  RAC  I  I  I, 1  =  1, NR  ) 

RIT  2,2,  (RACQI (),  1=1, NT) 

RIT  2,2,  (RHOm,I  =  l,NT> 

CALL  LRID  <662,  ARC  .  RAO  I 
C  READ  IN  AND  STORE  DATA  ONE  FILE  PER  TEMPERATURE 

72  NFILL=l 

cu  ieo  N  = 1 , N T 
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STUDY  OF  KADI  AT  ION  HEAT  FLUX  FROM  HIGH  PRESSURE  AIR  ARCS  AT  LCW  F09/11/64 
IF  (N-l)  73C.  73,  83 

C  READ  BCCREC  FDR  FIRST  FILE  ON  AS 
C  NW BCD  IS  SET  TO  2ERO  TU  READ  ENTIRE  RECORD 

73  ERASE  NWBCC 

IF  (RREADI645,BCCREC,NWBC0) 175,  74,75 

74  CALL  ERROR  I 29H0BCCREC  SHOULD  START  NEW  FILE) 

75  NPP=3 

76  CO  77  K=1,NPP 
ERASE  NWB  IN 

IF  (RREAO<66l,R6C,NWBIN ) ) 77, 78,77 

77  CONTINUE 
GO  TO  85 

7e  WOT  10.79,  NFILE, K 

75  FORMATI 13HCFILE  NO.  IS  I2,2XUHKEC  NC.  IS  12) 

CALL  ERROR  I 18HE0F  WHILE  READ  REC ) 

83  CALL  CLOSE  (661) 

NFILE  -  NFILE  *1 

NPP*2 

GO  TC  76 

85  IF!  ABSFIT£MP(N)-RECI 1))-1.E-6)90,9C,83 
9C  K=NPP+1 

ERASE  NWBIN 

IF  ( RREAD I  661 ,REC , NWS  I N ) 194,78,94 
95  IF  (ABSFIRhO(N)+2« )  -  l  • £-6 )  110,110,95 

95  IF  (RH0(N)*2. )96, 110,97 

96  PCNT  sRHOINJ+3. 

1  PST  =  l 

GO  TC  LOB 

97  l FI AHSFC  RHC IN) l-l.E-6) 14C, 14G.98 

98  IF  (RHCIN) 1100, 140, IC5 

IOC  IF! A8SFI RHCIN )*1. l-l.E-61 125, 125,101 

101  JFIRFCIN)*!. ) 102,125,1C3 

102  PCNT  =  RHOIN)  +2. 

IPST  =  2 

GO  TO  108 

103  FCNT  =  RHOIN )  +  1. 

IPST  *  3 

GO  TC  106 

105  IF  IABSFIRHOlNI-1. I— 1 ,E-6 ) 15 5, 155 , 106 

106  IF  IRHC(NI-1.)107,155,17C 

107  PCNT  =  RHOIN) 

IPST  =  4 

108  CU  109  I  =1, Nw 

PROP  IN,  I  )  =  EXPFKLOGFIREC  I  IPST+D/REC  I  IPST)  )  )  «PCNT*LOGF  I  REC  I  IPS 

IT)  )  ) 

105  IPST  =  IPST+5 
CO  TC  1  HO 
11C  IPST  *  2 
GO  TC  160 
L25  IPST  =  3 
GO  TC  160 
14C  IPST  =  4 
CO  TC  16C 
155  IPST  =  5 
160  CO  166  1*1. NW 

PROPIN, I)  =RECt  IPST ) 
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165  I  PST  =  IPST  +5 
GO  TC  180 

17C  PCNT  =HHCIN)-1. 

IPST  =5 
CO  175  1=1. NM 

PROPIN, II  =  EXP FI (UOGFIRECI IPSTI/RECIIPST-l) 1 1 *PCNT  + 

1L0GF1RECI IPST)) ) 

175  IPST  =  IPST  *5 
18C  CONTINUE 

CALL  CLOSE  1661) 

«0ALPH=DALPHA*1.7453292E-2 

NH=NDItf-l 

NANGLE  =NKAVE*INALPHA-1)+1 
ALPHA  I  11  =  0. 

SINAI I)=C. 

CO  190  NA=2 , NALPHA 
ALPHA ( NA )  =  ALPHA  I  NA-  l )«-R0ALPH 
SINA(NA)  =  SINFI ALPHA (NA)  I 
19C  CS I NA I NA )=SINAINA)-SINA INA-1 ) 

DO  210  NG-I.nTR 

CO  210  NA=»2 .  NALPHA 

EXKGSNI NC.NA )=RADINQ )/S IN A  I N A ) 

2 1C  EXSNRCINO.NA )  =  1./EXRCSNINC,NA) 

ANGLE! l)=0. 

ANGACO=RCALPHy H AvE 
CO  22C  NAVE=2,NANGLE 
ANGLE  (NAVE)=ANGL61NAVE-11 «ANGAOC 
220  EnSNA(NAvE)=SINFIANGLEINAVEI I 
I'.O  2  30  NQ=l,NT 
CO  230  NA VE=  2 • N ANGL E 
ENRGSN! NO, NAVE )=RADO(NG)/ENSNA! NAVE J 
230  ENSNRCINO.NAVEIsl./ENROSNINO.NAvE ) 

CG  2*0  N= 1, NT R 

24C  ABEKISI N)=-l .91  I L.-l RACING  1) /RAOINI )**2/2. )-(ARSINIRADIN+l)/RACIN 
1  I ) +SINFI 2.-ARSINIRA0IN+1) /RAO(N) ) 1/2.)/ 3. 1415926) 

CU  2*1  N=1,NT 

TEMPKH(N)  =  (TEMPINI*»4)«  5.6686E- 12 
RDELTAtNI=RAOIN  l/OCLTA 
CR<M  =  R4DIM-RA0IN«l) 

C2T«N)=t.438/TEMPIN) 

241  CV(N>=AINC*C2T(NI 
ERASE  C T A 
CC  7C0  1=1. NM 
ViNA  =  WNA-A1NC 
MM  I  I  =MNA 

ERASE  AMIN, A MAX. OPT .SUM.CRK 
CU  245  N-=  1 , N T 

245  SUM=SUM*OR(N}*PROP(N, 1) 

CK*2.*5UM 
EXT  RPN  =  1  . 

EXTRHX-1. 

CU  650  N= 1. NT 

EKOKI l)=OR(N)*PKOP(N.l) 

CRK  =  CRK  ♦EKCRI l ) 

V  =  Whl I ) *C  2  T I N ) 

IF  (V-. 278)650, 253,250 
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25C  IF  IV-L3. 3251,252, 65C 

251  IF  I V- I • 1253,253,252 

252  ETAIN, [1  =  . 15399 *DV I N 3  * ( V«* 3 ) /( EXPF I V )- 1 . 3 
GO  TC  255 

252  V2=V*V 
V  3=V2*V 
V4=V3«V 
V6=V4*V2 

X  =  V  2-V3/2 • *V47 12.-V6/720. 

2=(QVIM**21*C1./12.-V  /8.+V3/2A.-VA/5  76.  I 
ETAINfl  )  =  .  15399*DV(N)*(X*Z  I 

255  IF  IN-NT  >  256,650, 73C 

256  ERASE  SKOR 

CO  275  NA=2,NALPt-A 
LO  27C  NQ=l,N 

R(NU)=CXKgSNINO,NAl*ISlNF( AL PHA I N A 3-ARSIN I  RAD  IN* 1 1  * 
lExSNRCINO, NA)  )]  ) 

IF(  N  -1)  730,265,257 

257  IF  (NG-L)  730,270, 26C 

26C  AKCRINC-1 )» {RtNa-1 J*PROP(NQ-l, I 3 
IF  INO-  N  )270, 265,730 
265  flKCKI  A  l=R(  N  ) #PROP(  N  ,1) 

EKCK  I  NA  )  *  AKDK  [  N  J 

270  CONTINUE  / 

CO  275  •'J0=  l , N 

275  SKDR(NAI=SKDRINA)*AKCR<NCI 
SKUKI  l ) -CHK 
CO  305  NA  =  l  * NALPHA 
I  FI CXDK I NA }- 1 . I  295, 255, 3C0 
295  FKK(NA1=EXPF(-1.206*EKCR(NA]  I 
CU  TC  305 

JCC  IF  (EKCRI NAI-16. I  301, 3C1«302 

301  FKH (  N  A I  =  ( EXPFt— 1.2G6*EKDRlNA) I J  *  I  1 . 0774-. 09539* EKDR [NA I ♦ 

1-0555  77 • ( EKDRI NA I  I **2-A .29 1A3E-3*! EKOR I NA ) 3 **3*2.082602E-4* 

2IEKURINAI ) **4) 

-  GO  TC  305 
30  2  FK*<NA)=0. 

305  CO.nTINUC 

ERASE  CPSMN* EPSMX 
CO  310  NA=2, NALPHA 
EPSMM=EPSKN+FKK(NA J*CSINAINA) 

3 1C  GPSPX=£PSKX*FKR(NA-I 1*CS1NA( NA) 

EGUr=l«— EXPFIABEHISIM*  PROPIN, I)«RAOINJ  3 
E  IN= l.-I EPSHN+CPSMX) 72. 
fcBD-ETAIN,  I3*TEMPRNtN) 

VIIN  =  WBB«EIN*ROELTAIN*ll 
wUUT=ViCe»ECUT*RDELT  AIN) 

WOUTPN- WGUT  *EXTRMN 
WCUTPX=WCUT*EXTRMX 
CG  340  NA=1, NALPHA 
IF  ISKCRINA3-1. 3320, 32C, 325 
320  FKR<NA)=EKPF(-1.206*SKC8INA) 3 
GO  TC  34C 

325  IF  ISKCRI NA3-16. )  330,330,335 

3  30  FKR I NA) = I  EXP  F I  - L . 2C6*SK0R ( NA I )) *  1  1 . 07 74-.09539*SKDR I NA ) ♦ 

1.055577*1 SKDRINA) )**2-A . 29 143E- 3* < SKDR I NA } )**3+2.082602E-4* 
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2(SKOR(NA))««4) 

GO  TC  340 
335  FKFMNAJ  =0. 

34C  CONTINUE 

ERASE  6XTRHN,  6XTRMX 
CO  350  NA=2,NALPHA 
£XTRMN=EXTRMN*FKRtNA )*C SI  NAINA) 

35C  EX  T RMX=EXTRMX  *FKR ( NA- 1 J  *0 S IN A< N A ) 

ERASE  TRYKCR 

NT£ST=NT-N-1 

CO  6C0  NAVE=2,NANGLE 

ERASE  K  I 

ERASE  NTCUCH 

CO  530  NG= 1, NT 

S  1  GMA=RAD(  N+-1 1  *ENSNROI  NQ«  NAVE  1 
IF  IStGMA-l. 1519,519,518 
518  MOUCH  -NO 
CO  TC  530 
515  B=ARSIMSIGMA] 

IF  (NQ-NTEST I  520,520,523 
52C  C=ANGLE(NAVE)-B+3. 14I5S26 
I F I D- 1 . 5707963)  522,522,521 

521  C  =  3.  I415926-C 

522  HIINC,21=ENRCSN<NQ,NAVEI*SINF(0) 

522  C- ANGLE (NAVE  I *B 

I F ( C- 1*5707963)  525,525,524 

524  C=J. IM5926-C 

525  R1 (NC.1I*ENRCSN<NQ,NAVE>»SINF<C> 

530  CUNTINUE 

NFAC-NTOUCH  *-1 
NJ  =  Q 

IF(NFAC-NTEST) 545, 545,565 
545  CO  560  NG=NF AC, NTEST 
NJ=NJ*l 
J=NT-NC 

56C  AKDR(NJ)  =  (RHN0,2  )-Kl  CNQ*1.2  ))*PROP(  J,I  > 

565  N J=NJ ♦ l 

J-NT-NTOUCH 

AKOK(NJI  =  CRUNFAC,  1I-R  I  ( N  F  AC  ,  2  )J  «  PROP  I  J  ,  I  I 
IF  l NFAC-NTR  )  580,580,595 
58C  CO  590  NQ-NFACpNTR 
NJ-NJ+1 
J=NI-NC 

59C  AKDR(NJ)  =  (RI(NQ<-l,i)-RHK0,M»«PR0P1  J,(  ) 

595  NUMDR=fvJ 

CU  6C0  N J= 1 , NUMOR 

&0C  TRYKDKINAVE)-TRYKDR(NAVE) ♦AKDRINJ I 
606  TRYKCR ( 1 )=CK-ORK 
MAL-l 
NAL  =  1 

CO  635  LAL=l,NQIV 

NAL=NAL*NRAVE 

TEHPL  =  TRYXCR(NAU 

JAL=NAL-1 

CO  63C  Il=MAL,JAL 

IP l = 1 1  + 1 
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STUDY  OF  RADIATION  HEAT  FLUX  FftGH  HIGH  PRESSURE  AIR  ARCS  AT  LCM  F09/11/64 
CO  630  I2=IP1,NAL 

1 F I TRYKUR ( It l-TRYKDRI 12  1  )6 30 , 6 30 , 62 T 
627  TEMPS=TRYKCR(  11) 

TRYKCRl I  1 ) =TRYKOR( 12) 

TRYKCRl r2)=TEHPS 
63C  CONTINUE 

SKCRI LAL )=TRYKDRIMAL) 

IAL  =L AL  *NC  IV 

SKCRI IAL)=TRYKDRINAL) 

PAL=K AL+NRAV  E 

635  TRYKCRl MAL )=  TEMPL 
CO  64C  NA=l, IAL 

IF  I$KDRlNA)-l. 1636,636,637 

636  FKRINA)=EXPFl-1.206*SKCR(NA) J 
GO  TO  640 

637  IF  ( $KCR ( NA ) -  16 ♦ )  638,638,639 

6  38  FKR I N A 1  = I EXPF l- l  .206  *SKOR 1 NA ) ) ) * 1  1 . 0 774-. 09539*SKDft I NA 1  * 

1.055577*1  SKORINA) )  «»2-4 .29  14  3E-  3*1  SKCRI  NA)  )**3+2 . O026O2E-4* 

2 1 SKOR I NA ) ) »  *  4 ) 

GO  TC  640 
635  FKR(ISA)=0. 

64C  CONTINUE 

ERASE  ENTRPN, ENTRMX 
CO  645  NA=2, N ALPHA 
NAH=NA*NH 

E  NTRFN^ENTR^n  tFKR I  N  AH  )*CS1NA(NAI 
645  ENTRMX=£NTRMX4FKRl N A- 1 1 «CS IN  A I NA ) 
kINMN  =  ViIN*£NTRMN 
h INMX-KIN*CNTRMX 
AMIN=AHIN*WINMN4W0UTPN 
AMAX=ARAX4HINNX*HuUTHX 
CPT=CPT  4W  lN  +  kiOUT 
65C  CONTINUE 

CPTCCKI  I  l  =  TEMPRW<  NT  1 «RCELTAI NT) »ETA{NT, I ) *  I  1 .-E XPF t- 1 . 9*EKDR( 1 > ) 1 
COrtMM I  1 =UP  rCOR l I ) *EXTRMN 
CQRMX I  I  I  =OPTCO«I  I)*EXTKMX 
CHIN! I ) =AM I N+CQRMNl I > 

C  M Ax l l 1  =  AMAX  4C0RMX l I ) 

7CC  CCJPTI  I)=OPT+OPTCOR<  l) 

IF  I  1  CP Tl)  7 05, 703, 705 

703  CO  704  1=1, NR 
ERASE  Y 

Y  I  1 )  =ViN  l  l  1 
Yl 2)  =  C,MAX( I ] 

Y I  3 ) =0P  I N ( I  I 
Y 1 4 ) =  OGPT (  I  ) 

704  CALL  RITE  1662, Y, 4) 

CALL  CLOSE  1662) 

705  ViOT  1C,  l 
hOT  10,8 
WCT1C. 7 

CC  TG7  1=1,  NT 

707  VOT  10,  9,  TEMPI  l ) ,  RHO  <11,  RADII) 

V»DT  1C.1 
v»or  io,4 
CO  710  1=1, NW 


113 


AEDC-TR-65-1 l 


STUDY  OF  RADIATION  HEAT  FLUX  FRQH  HIGH  PRESSURE  AIR  ARCS  AT  LGM  F09/11/6' 

TIC  hOT  10,5#  MNIll.  OMAXtl)  ,  QMIN(l),  OCPTII) 

ERASE  QSUHMN,QSUNHX, QOPTSM, SMCRNN , SHCRMX»SHOPTC 
CO  711  1*1, NH 
CSUKHN=QSUFHN+QM1N< I ) 

CSUN«X=QSUMMX+QHAX( I ) 

COPTSM  «QOPTSM*QOPT( I  ) 

SHCRMN=SMCRKN+CORHN ( I ) 

SMCRHX=SMCRHX+CORHX I  I ) 

711  S HOPTC* SHOP TC+OPT COR ( I ) 

WOT  10, 13,SMCRNX» SHCRMM , SHOP  TC 
ROT  10,6,  CSUMMX,QSUMHK,QOPTSH 
ERASE  ETASH 
DO  722  N=  1 » NT 
CO  722  1*1, NW 

722  ETASK(M)*ETA(N,  l  H-ETASHIM 
ROT  10,10 
CO  725  1=1, NT 

725  ROT  10,11,  TEMPI II.ETASMI I) 

726  NCASE=NCASE-1 

IF  INCASE  I “55,740,72? 

727  CONTINUE 

F  P2E  661, , 32384 

GO  TC  71 

730  CALL  ERROR  I 10H INDEX  OUT  OF  LIHIf) 

74C  CONTINUE 

H  HIE  661 , , 32384 

H  HZE  662,, 32384 

CALL  EXIT 

END (1,1,0, 0,0, 1,1,1,C,1,C,0,C,C,C) 
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COMMON  ABEMIS(22),AINC,AKDR<  45  J  .ALPHA  <  10]  .ANGLE  (20) 

COMMON  BCDREC ( 25 ) 

COUPON  BUFF 2(  52 5 ) , BUFFT C 36 1 , CORMN ( 1 21 > ,CORMX( 121 > , C2T ( 23) 

COMMON  DALPHA.DELFA.CPHl 

CUHMCiN  DR  1231  ,DSINA<  10)  .EKDft  ( 101  .ENSNA  (20) 

COMMON  E.NSNRQ(2  3,20)  ,ENRQSN(  23, 20  )  , ETA  ( 23 , 121 1 
COMMON  E  T  ASM { 23 ) 

COMMON  EXSNRO( 22. 2C ) .EXRQSNl 22, 20 ) TFKR ( 20) ,FMT( 12) , t T (23) 

COMMON  (GPT1,  NALPHA ,ND 1 V 
COMMON  NRAVE, NTR, NT, NW,OPTCOR( 121) 

COMMON  PRCP (  23,  121 I.QMAXi  121  ),QMIN(  121) 

CUHHCN  QUPT 4 121), RA DO (23) ,RAD( 23 ) , RAVE , RDALPH,RDELTA ( 23 1 ,REC 1 700 ) 
COMMON  KHCH23) 

COMMON  RI(23,2) ,K( 221 ,S  INA ( 20) , SKDR (20 ) , TEMP (23) ,TEMPRH(23> 

COMMON  TRYKDRI20) t  feN A , UN ( 121 l,Y(4) 

1  FORMAT  I  IH 1  I 

2  FURHAT I  12A6J 

4  FORMAT  ( 1H03X6HWAVE  NO . , UX4HQMAX , 1 3X4HQMIN , 1 3X4HQ0PT ) 

5  FORMAK  IH  4  (  LPE  14.  7, 3X )  > 

6  FORMAT  I IH03X8HSQMAX  =  1PE 14. 7, 3X8HSGM IN  -  1PE14. 7t3X8HS0GPT  =  IPE1 
14./) 

7  FORMAT  (  IHC2XHHTEMPETATURE,  5X1 2H LOG  RHD/RHOO , 8X6HRAD1  US ) 

8  FORMAT! 1H030X19HT EMP ERATURE  PROFILE) 

5  FORMAT ( IN  31 IPE14.7, 3X) I 

1C  FURMAT  ( 1HC2X 1 1HTEMPERATURE, 8X6HETASUH) 

11  FORMAT  ( IH  2 ( 1PE 14.7,3X1) 

13  FORMAT  t 1H01X10HSQCORMX  =  IP  El 4. 7, 1X10HSQCCRMN  •  1PE14.7, 

1 IX 1 IHS CCCRCPT  x  IPE14.7I 
C  INITIALIZE  FOR  INPUT  OUTPUT 

5C  L  F I  IMT(642,643»8UFFT(86),  39, 47)  160,55,60 
55  CALL  DUMP 

C  INITIALIZE  FUR  ERROR  ROUTINE 
6C  IF  (ERTRP4 5)170,726, 726 
C  INITIALIZE  FOR  DUFFERS 

C  A6UF 1  IS  FOR  INPUT  ON  A5 

70  CALL  AEUFl 
RIT  2,2,  INCASE 

IF  { BUFFS ( 662, BUFF2( 525), 525,0,1,1,1)) 71,55, 71 
C  READ  IDENTIFICATION  CARC 

C  REAC  IN  CAROS  CN  A2 

71  RIPT2.2, I  FMT  ( I ) , 1=1, 121 
CALL  SPGhOR  I FMT  ) 

C  INPUT 

RIT  2 ,2  ,  CELT  A  ,VINA,  WNL ,  C  ALPHA  ,RAVE 

HIT  2,2,  NT,  IOPTI, 10PT2 ,NTR, NR , NALPHA »NRAVE »NDl V 

RIT  2.2, (TEMPI  I ) , 1  =  1, NTR) 

RIT  2,2,  ( R AO ( I ) , I = 1 ,NR ) 

RIT  2,2,  (RAOO(  I), I  =  1 , N  T ) 

RIT  2,2,  (RHOtl > , I  *  1  ,NTR ) 

ERASE  PRCP 
NFILfc  =0 
ERASE  NWBCC 

IF  (RRE AO I645.BC DRECtNW BCD)) 75,74,75 
74  CALL  ERROR  (25H  BCOREC  SHOULD  START  TAPE) 
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75  ERASE  NWBIN 

IF( RR EAD t 66 1,REC, NWBIN) 176,74,76 

76  CO  18C  N- I, NTR 

77  NFiLE  =NF ILE ♦ 1 
K  =  1 

ERASE  NWBIN 

IF(RREAD(661,REC,NWB1N) >83,78,63 
7e  WOT  10, 79 » NF ILE » K 

75  FORMAT l I3H0F ILE  NO.  IS  12, 2X 1 IHREC  NO.  IS  12) 

CALL  ERROR  ( I8KE0F  WHILE  READ  MEC) 

83  IFIABSFI TEHPINl-RECI 1) )-l.E-6)90»90»65 
85  CALL  CLOSE  (661) 

GO  TO  77 
9C  K  =  K*l 

ERASE  NWBIN 

IF  |RREADt661,REC,NHBIN) J94, 78,94 

94  NW=RECU) 

l F I N-NT  R 1192,  1190,  192 
1I9C  CO  191  1=1, NW 

191  UN(  I  )=RECI  1+21 
hN I NW  + 1 1  =  HNL 

192  K  =  K+ 1 
ERASE  NWS  IN 

IF  IftftE ADC 661, REC, NWBIN ) 1 194,78,194 
194  IF  I ABSFCRH0IN>*2. 1-1. £-6)110, 110, 95 

95  IF  (RHCtNl*2.)96, 110,97 

96  PCNT  -RH0(N)*3. 

1PST  *  1 

GO  TG  108 

97  IFCABSFCRHCCNI l-l.E-61140, 140.98 
96  IF  CRHOINDIOO, 140,105 

100  IFIABSFI RHCIN1+1.) -l.E-6) 125, 125, 101 

101  IFIRhOINl+l. 1102, 125, 1C3 
1G2  PCNT  =  RHOIN I  *2. 

I  PS  T  =  2 
GO  TO  LOB 

103  PCNT  =  RHOIN)  ♦  1. 

I  PST  =>  3 
GO  TO  100 

105  IF  (A8SFIRH0INJ-1. J-l.E-6) 155, 155,106 

106  IF  (RHCtN)-l. 1107, 155, 170 

107  PCNT  =  RHO(N) 

IPST  =  4 

10B  CO  109  I  = 1, NW 

PROPIN,  I)  =  EXPF(CLOGF(REC(IPSl*l)/RECIIPST)l  I  •PCNTHOGF  I  REC  I  IPS 
IT  )  >  ) 

109  IPST-  =  IPST*  5 
GO  TC  180 
11C  IPST  =  2 
GO  TO  160 
125  IPST  =  3 
GO  TG  160 
14C  IPST  *  4 
GO  TG  160 
155  IPST  ■  5 
160  CO  165  1=1, NW 
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PROPIN. II  *RECI  IPST I 
165  I  PST  »  IPST  *5 
GO  TO  LBO 

170  PCNT  -RHOI N ) -1 . 

IPST  *5 
CQ  175  I=1,NW 

PROPIN. I  I  •=  EXPF 1 1  LOG Ft  RFC  I IPST1/RECI IPST-1 ) I I*PCNT  *  . 

1LOGFIRECI  IPST))  I 
ITS  IPST  =  IPST  45 
I8C  CALL  CLOSE  1661) 

IF  1 1 OPT 2)184.193)184 

183  NTA-NTR+1 
NTB=NTR 
GO  TC  185 

184  NTA-NTR 
NT0-NTR-1 

185  NTP-NT-1 

RDALPHaOALPHA«l.  7453292E-2 
K'H-ND  I  V~1 

NANGLE  =NRAVE»INALPHA— 1 )+l 
ALPHA! 1>*0. 

SINAI  11=0. 

CO  190  NA  =  2. NALPHA 

ALPHA  I NA )- ALPHA  I N A- 1 )+RDALPH 

SINAINA]=SINFI ALPHA! NA ) ) 

19C  CSINA(NA)>SINAINA)-SINA(NA-1 ) 

GO  210  NQ*1,NTB 
CO  210  NA=2, NALPHA 
EXHUSNING.NA )»KAD( NO )/S INAINA ) 

2 1C  EXSNRGINQ,NA)=1./EXRQSNINQ,NA) 

ANGLE  1 1 > =0. 

ANGACD-RDALPH/RAVE 
DO  220  NAVE-2, NANGLE 
ANGLE  I  NAVE  I =ANGLE I  NAVE-1 1 4ANGADD 
220  ENSNA (NAVE)=S1NF( ANGLE! NAVE)  ) 

CO  230  NQ=1,NT 

CO  230  NAVE=2, NANGLE 

ENRQSN I NG . N  AVE ) =RADQ ( NG ) /ENSNA I  NAVE  > 

230  ENSNR0INQ«NAVE)=1./£NRQSNI NO, NAVE > 

CO  241  N- l.NTR 

A0EM I S ( N )*-l .91 ( l.-(«AD(N4l) /RADIN) l**2/2. )-l ARSINIRADIN+IJ /RAOIN 
1) )4SINFI2.«ARSINIRA0IN4l)/RADIN)) )/2.)/3. 1415926) 

TEMPRNI N)={TEMPIN) *»4| *5.66B6E-12 
RDELTA(N)*KAOIN l/OELTA 
CR(N>*RADIN>-RAD(N4l) 

241  C2T(NI=l.43B/(TEMP(N)*4.) 

ERASE  ETA 
CO  7C0  I-WNH 

ERASE  AKIN, AM AX .OPT , SUM.DRK 
CQ  245  N-l.NTR 
245  SUM=SUM4CRIN)*PR0PIN,ll 
CK-2. *SUM 
EXTRMN=1. 

EXTRHX=l. 

CO  650  N- 1 , NIA 

CV=IHN(  IH)-WNA)*2.*C2T|N) 
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EKOR'1 1  )ZDR(NI  •PROP  IN*  1 1 
CRK=ERK+EKCR l 1) 

V=C2TIN)*(2. »WN I  I ) ♦HNA+HN ( 1+ I ) I 
IF  tV-13. 1252,252. 65C 

252  ETAIN, I)*.15399*0V  *( V**3I /IEXPFIV1-1. ) 

255  IF  IN-NTA)256,650, 730 

256  ERASE  SKOR 

CU  275  NA*2. NALPHA 
CU  270  NC=1,  N 

R(NQ)*EXRQSNINQ,NA)*tS  INF  ( ALPHA  (KA)-ARSIN4RAD(N+M« 
lEXSNRQtNQ.NA) ) ) ) 

IF!  N  -1)  730,265,257 

257  IF  (NO-II  730.270.26C 

260  AKOKt  NO-  l)=CRtNQ-U-RtNQ)  )  »PROP  t NQ- 1 ,  1 1 
IF  tNO-  N  ) 270, 265, 730 
265  AKDRt  N  )~R(  N  )*PROP<  N  ,1) 

E  KDR I NA ) ~ AKDR (  N  ) 

27C  CONTINUE 

CO  275  NO=I, N 

275  SKOR|NA)=SKOR(NA) +AKDR  ( NO  I 
28C  SKDR ( 1 )  =  ORK 

CO  305  NA  »  L , NALPHA 
I F  t  EKORt  NAI - 1 . ) 295 , 295, 3C0 
295  FKR  CNA I  =  EXPFC -1 . 206»EKCR I NA) ) 

GO  TO  305 

30G  IF  ( EKCrt (NA)-16.)3Cl,3CIf 302 

301  FKR I NA J  =  ( EXPFt  —  L.206»EKDR { NA )  )  )*t  1 . 0774-. 09 539» EKDR I NA ) ♦ 

1.055577*(EK0RtNAI J **2-4 .29143E-3* ( EKDR t NA M *• 3*2. 082602E-4 » 

2 ( EKDR (NA ) )**4 1 

GO  TC  305 

302  FKRt NA J  =  0. 

305  CONTINUE 

ERASE  EPSHNtEPSMX 
CO  310  NA*2,NALPHA 
EPSHN=EPSMN*-FKR(NA  i *cs inaina  > 

310  E  PSMX=EPSMX  *-FKR (NA-1)*C5INA(NA) 

EOUT»l.-EXPF(ABEMlSIN I* PROP! N, 1 1 *RAO(N I  ) 

EIN-l.-(EPSMN*EPSMXI/2. 

»(BB=ETAtN,I  l*TEHPRWtN) 
to IN=HBB*E IN*ROEL  TA { N  + 1 ) 

HOUT=WBB»ECUT •RDELTA(N) 

KOUTMN*WOUT*EXTRMN 
RaUTMX-WaUT*Exr«MX 
CO  340  na*i,nalpfa 
IF  (SKGRINAl-l. I32G, 120,325 
320  FKK(NA)  =  EXPFt-1.206*SKCR(NA)  ) 

GO  TC  340 

325  IF  { SKCH  t  NA ) - 16 . ]  330,330,335 

330  FKRtNAIM EXPFI - 1 . 206»SKOR t NA ) J )•< 1 . 0774-. 09539*SKDR ( NA ) ♦ 
1.055577*(SKDR(NA)1 **2-4 .2 9 14 3E- 3* < SKCR t NA » ) **3+2.082602E-4* 

2(SKDR{NA))**4) 

GO  TC  340 
335  FKRt  A A  J -  0 . 

34C  CONTINUE 

ERASE  EXTRMN.EXTRMX 
CO  350  NA=2, NALPHA 
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EX  TRMN=E  XTRMN+FKR INA )«DS I NAI  NA ) 

350  EXTRMX=EXTRMX+FKRINA-1) *CSINA|NA) 

ERASE  TRYKDR 
NTEST=NT-N-1  . 

CO  600  NAVE=2,NANGLE 
ERASE  Rt 
ERASE  NTOUCH 
CD  530  NQ= It  NT 

S 1 GHA=R ADI N+1)*ENSNRQINQ, NAVE) 

IF  I  S  IGMA- l .  ) 519*519,518 
518  NTOUCH  *NQ 
GO  TG  530 

515  13  =  ARS  [MS  IGMA  ) 

IF  (NQ-NTEST 1520,520,523 

520  C=  ANGLE r  NAVE  J-BO.  1415926 
IFI 0-1. 5707963)  522,522,521 

521  C=3. 1415926-D 

522  R 1 1  NQ » 2 ) =£NRGSNi NQ, NAVE )«SINFID1 

523  C= ANGLE C  NAVE  3  +B 

IFIC-l. 5707963)  525,525,524 

524  C=3. 1415926-C 

525  RI INC, 1) =ENRQSNINQ,NAVE l*SINFIC) 

530  CONTINUE 

N  FAC=N  FOUCh* 1 
N  J  =  0 

IF  INF  AC-NTEST ) 545,54  5,565 
545  CO  56C  NC=NFAC, NTEST 
N J  =NJ  +  1 
J=NT-NC 

560  AKOR <NJ)=lRItNQ*2)-RIINQ+l»2))*PROPIJ,II 
565  NJ^NJ+1 

J=NT-NTOUCH 

AKDR(NJ)=IRt INF AC, 1  )-R  I INF AC, 21 ) *PRQP I J, I ) 

IF  I NF AC-NTP 1530,580,555 
58C  CO  590  NQ=NFAC,NTP 
NJ=NJ+1 
J=NT-NC 

59C  AKDR I N J )  =  I R 1 1 NQ ♦ I , I ) -R I ( NQ , 1 > ) *PR0P ( J,  I ) 

595  NUMDR=NJ 

CO  600  N J=  1 , NUMDR 

60C  TRYKCRINAVE)=TRYKDRINAVEI+AKORlNJ) 

606  TRYKCRI I )*CK-0RK 
KAL=1 
AAL=  l 

CO  635  LAL=1 , NDI V 

NAL=NAL*NRAVE 

TEMPL=TRYKDRI NAL I 

JAL=NAL-1 

CO  630  1 1=FAL«  J AL 

IP1=I  1U 

CO  63C  I  2= I P 1, NAL 

IFI TRYKDR I  1 1 l-TRYKCK I  1 2 l > 630 , 630 , 62 7 
627  TEMPS=  TRYKCR (II) 

TRYKCRI 1 1 )  =  TRYK0R I  12) 

TRYKCRI I2)=T6HPS 
630  CONTINUE 
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SK0K(LAL)=TR7XDR|NAL> 

IAL  =  L  AL ♦NCI V 

SKORI IAL >=TRyKOR(NAL ) 

MAL=HAL*NRAVE 

635  TRYKDR(MAL)=TEMPL 
CO  640  NA=l, IAL 

IF  (SKCR1NA 1-1. >636,636,637 

636  FKR(NA)=EXPF(— i«206*SKCR(NA) } 

GO  TC  640 

637  IF  ( SKCR( NA ) -16. I  638,638,639 

636  FKR ( NA 1= ( EXPF(-1»206*SKDR( NA ) ) I  * ( 1 . 0T74-. 09539* SKDR {  NA)  ♦ 

1. 055577 •(SKDR(NA) >**2-4.29 143E-3*(SK0RINA ) )**3+2.082602E-4* 

2( SKDRINA I ) **4 1 
GO  TC  640 
635  FK«(NA)=0. 

64C  CONTINUE 

ERASE  EN 1KKN. ENTRHX 
CD  645  NA=2, NALPHA 
IS  AH=NA*NH 

eNTRf^ENTRMNtFKRl  NAH  1*DS1NA(NA) 

645  ENTRMX=ENTRMX*FKRINA-U •OSlNAINA) 

WINMN*WIN«ENTRMN 
U INMX=W  CN»ENTRMX 
AMIN'Af'  IN+WINMN*WOUTMN 
AMAX= AMAX+K INMX+WOUTWX 
Cpr=cpT+wiN*nour 
650  CONTINUE 
WNA=WNt I ) 

IF  tICPT2>67O,680,67C 

670  CPTCGRI II*TEMPKW(NTR>*RDELTAINTR)«ETA(NTR,I 1 
1*(  l.-EXPF[-l.9*EKCR<  111! 

CGRMNI  lUOPTCORI  I  >  *EXTRMN 
COKMXm=OPTCOR(  I)*EXTRMX 
680  CMINI  H  =  AMIN*CORHN(  l ) 

CMAX(  1 )*AMAX*CDRHX(  1  ) 

70C  COPTI l  )*OPT*DPTCOR< 1 ) 

IF  ( IOPTl)705,703,705 

703  CO  704  ( =  1 ,  NW 
ERASE  V 

7{ i)=WN(  I  ) 

7(2) =QMAX ( I ) 

7  (  3>  =  0FMNI  I  ) 

7  I  4 )  =  OCPT ( I > 

704  CALL  RITE  1662,7,41 
CALL  CLOSE  (6621 

705  SOT  10,1 
WOT  10, B 
UUH0,7 

CO  707  I  - 1 ■ N I R 

707  SOT  10.  9,  T EMP (It,  RHCU),  RADII) 

SOT  10,1 
SOT  10,4 
CO  710  1=1, NW 

7 1C  SUT  10,5,  SN(l),  QMAX ( I  ) ,  OMIN(I),  QQPT(l) 

ERASE  CSUMMN,QSUHMX,CDPTSH,SMCRMN, SHCRMX , SMOP TC 
IF  U0PT2J1710, 712.171C 
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1 7 1C  CO  711  1=1, NW 

SMCRHN=SNCRMN»CORHNI  II 
SMCRPX=SMCRMX*CORHX(  (] 

711  SMOPTC=SMUPTC*OPTCOR(I) 

712  CU  715  1=1, NW 
CSUMNN=Q$UMMN»QMIN( I  ) 

CSUMMX  =  QSUMMX-KJMAXI  I  ) 

715  COPTS M=QOPTSM+QOPT l  I  1 

hOT  1C, l 3, SMCRMX, SMCRHN , SMOPTC 
WOT  10,6.  GSUMMX , QSUMMN , COPT  SH 
ERASE  ETASK 
CO  7 22  N=1»NTR 
CO  722  1=1, HU 

722  ETASMCN»=ETA(N. Il+ETASKINJ 
hGT  10,10 
CO  725  1=1, NTR 

725  hOT  10,11,  TEMPI  n.ETASMU) 

726  NCASE-NCASE-1 

IF  INCASE)  55,750,727 

727  CONTINUE 

H  H4E  661, ,32385 

GO  TC  71 

73C  CALL  ERROR  I18HIN0EX  OUT  OF  LIMIT) 

75C  CONTINUE 

H  MiE  661, ,32385 

H  M2E  66?,, 32385 

CALL  EXIT 

END (IflrQtGtOt l  * li 1 f  Ct ItCpOiCtOpC) 
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APPENDIX  C 

CALCULATION  OF  ANNULUS  RADII 


C.l  ANALYSIS 


If  one  assumes  that  an  individual  annulus  is  optically  thin,  then 
there  exists  a  mean  value  J  of  total  spectral  radiance  (watts  /cm^  - 
ste radian)  such  that 


4  it  J 


,2  2  . 

"(rj  -rj«\ 

•'l 

Zn  r  (4tt  J)  dr 


j+1 


(C-l) 


where  Eq,  (C-l)  is  written  for  a  unit  length  of  arc  column.  Further 
assuming  that  J  is  proportional  to  and  knowing  T  as  a  function  of  r, 
Eq.  (C-l)  can  be  Rewritten 


T 

J 


(C-2) 


For  measured  profiles,  the  function  T(r)  was  approximated  as  a 
polynomial  using  a  standard  "least  squares"  curve  fit  computer  pro¬ 
gram  (Appendix  D).  Starting  at  the  outer  boundary  of  the  arc  column, 
successive  values  of  rj+j  were  determined  such  that  Tj  was  a  tem-Q 
perature  at  which  absorption  data  were  tabulated  (multiples  of  1000  K). 


Since  appears  both  explicitly  and  as  an  integration  limit, 

the  following  computerized  iteration  scheme  was  used.  The  outside 
and  inside  radii  of  the  j**1  annulus  are  r  ■  and  r^  respectively  and  the 
radius  at  which  T  =  Tj  is  r  (Tj)  (see  sketch).  The  radii  r  (Tj)  and 
r(Tj+j)  were  chosen  as  trial  values  of  rj^.j  and  the  weighted  average 
temperature  computed  from  Eq.  (  C-2).  Additional  successive  trial 
values  of  rj+^  were  obtained  by  linear  interpolation  from  the  two 
preceding  trials.  Typically  four  trials  were  needed  to  match  the 
weighted  average  temperature  to  the  tabulated  T-  within  one  degree 
Kelvin. 
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C.2  DESCRIPTION  OF  ANALYTIC  TEMPERATURE 
PROFILE  COMPUTER  PROGRAM 

A  FORTRAN  program  called  ANALYTIC  TEMPERATURE 
PROFILE  was  written  for  the  IBM  1620  Digital  Computer  and  used  to 
calculate  the  successive  values  of  r.  .  satisfying  Eq.  (C.2).  It  re¬ 
quires  as  input  data  the  constants  of  “me  polynomial  describing  the 
profile  {see  Appendix  D)  and  a  list  of  the  Tj*. 

Profiles  derived  from  experimental  data  generally  were 
incomplete.  That  is,  the  data  extended  neither  to  the  centerline  nor 
to  the  outer  boundary.  This  program  therefore  automatically  com¬ 
pleted  the  profiles  with  parabolas  which  joined  the  curve  fit 
polynomial  smoothly  and  also  satisfied  the  conditions  of  zero  tem- 


*  o  o 

From  10,  000  K  up,  increments  of  1000  K  in  T:  were  used  to  make 

maximum  use  of  available  data.  Below  10,  000°  K  increments  were 
usually  2000  K  because  in  this  temperature  range  the  temperature 
gradient  was  generally  very  steep  and  the  1000  K  increment  would 
result  in  annuli  unnecessarily  thin  compared  to  the  higher  tem¬ 
perature  ones. 
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perature  gradient  at  the  centerline  and  a  given  outer  boundary  tem¬ 
perature.  Next,  a  calculation  of  the  values  of  r{Tj)  was  necessary. 
Since  the  curve  fit  polynomial  was  usually  at  least  3rd  order  the 
Newton -Raphson  iteration  method  was  used. 

Arbitrary  profiles  of  the  form  of  Eq.  3.  3  were  complete  without 
end  parabolas  and  could  be  solved  exactly  for  the  r(Tj).  A  simplified 
version  of  the  Analytic  Temperature  Profile  program  was  used  with 
these  profiles.  The  statements  in  the  program  listing,  Section  C,  5, 
which  are  bracketed  were  deleted  and  four  statements  were  added  to 
compute  the  r(T.), 
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C.3  DEFINITION  OF  SYMBOLS  FOR  COMPUTER  PROGRAM 


A(N) 

B,  B2 

C,  Cl 
COEF 

DELT 
DRAD 
E,  El,  E2 
F(I) 

FINTG 

FPOW 

I 

UMAX 

K 

L 


Constanta  associated  with  curve  fit  polynomial 

Constants  associated  with  inner  parabola 

Constants  associated  with  interpolation  of  RLUM 

Coefficient  used  in  calculation  of  R(N)  by  Newton- 
Raphson  Method 

Number  of  intervals  for  numerical  integration 
Radius  increment 

Constants  associated  with  outer  parabola 

Function  of  radius  F(I)  =  jr{R)~]  ^  x  R 
r 

]  FdR 

After  some  calculation  steps  FPOW  =  T(R) 

Index  for  integration 
Maximum  vaLue  of  I 
Index 

Index  for  computing  R(N)  and  for  computing 
.fRJLLIM^L),  ’IAVG(L)] 


N 

NCST 

NORD 

NT 

PRNUM 

R(N) 

RAD 


Index  for  temperature 

Number  of  constants  in  curve  fit  polynomial 
Order  of  curve  fit  polynomial 
Number  of  temperatures 
Profile  number 

Radius  associated  with  temperature 
Radius 
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RB 

RD,  RE 

RDIF 

RJ(N) 

RPOW 

RLLIM(L) 

RTOL 

RUL-IM 

SUM1,  SUM2, 

T(N) 

TAVG(L) 

TD 

TDIF 

TDIFA 

TDP 

TE 

TEP 

TNR 

TP  NR 

TTOL 


Intermediate  variable  in  computing  R{N) 

Radii  at  which  inner  and  outer  parabolas,  respectively, 
join  the  curve  fit  polynomial 

Error  in  Newton-Raphaon  calculation  of  R(N) 

Inside  radius  of  annulus  having  a  weighted  average 
temperature  T{N} 

Intermediate  variable  in  computing  R{N) 

Lower  limit  radius  for  TAVG(L)  integration 

Tolerance  on  Newton  Raphson  calculation  of  R(N} 

Upper  limit  radius  for  TAVG(L)  integration 

SUM3  Intermediate  variables  in  TAVG  integration 

Nth  temperature 

Weighted  average  temperature 

Temperature  at  R  =  RD 

Difference  between  TAVG(3)  and  T|N) 

Absolute  value  of  TDIF 

Slope  of  temperature  profile  at  R  =  RD 

Slope  of  temperature  at  R  =  RE 

Slope  of  temperature  profile  at  R  -  RE 

Temperature  as  used  in  Newton-Raphson  calculation 

Slope  of  temperature  profile  as  used  in  Newton- 
Raphson  calculation  of  R(N) 

Tolerance  on  TDIF 
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Read  Special  Data: 
RD,  RE,  N'T,  NCRD, 
PRNUM,  DEJLT, 
IMAX,  RTOL,  TTOL 


Read  Temperatures 


ompute:  TD,  TDP, 
TE,  TEP 


Compute  radius  R{N) 
associated  with  each  T(N): 

1. 

T(N)iTE:  Compute  R(N) 
directly  from  11 E"  parabola 

2. 

TE  <  T(N)<  TD  solve  for 

R(N)  by  Newton-Raphson 
method 

3. 

TD^T(N)  compute  R(N) 
directly  from  "B"  parabola 

a p  6  T{N) 


itle,  special 
data.  TD,  TDP,  TE, 
TEP,  and  coefficients 
for  parabolas 


To  190 


TAVG(L) 


KRULIM^  -  (RLLIM(L)^ 


NOTE;  T(R)  is  given  by  the  curve  fit  polynomial  or  by  one  of  the  end 
point  parabolas,  depending  on  R 


F rom  229 


Yes  <RLLIM(3)  = 
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Yesc'ff AVG(3)  -  T 


i-T® 


TAVG(3}<  T(N>  L.-  No 


Set  RJ(N>  =  0 
for  current  value  of 
N  and  for  all  higher 
N  through  N  =  NT-1 


Yea  <TClLIM(3)^R(N  +  lj 


TAVG(2)  — *■  TAVG(l) 
RLLIM(2)  — >RLLIM|1) 

TAVG(3) - -  TAVG(2) 

RLLIMI3)— *  RLLIM(3) 


To  240 


RJ(N)  =  0 
N  =  N  +  1 


RJ(N)  =  RLLIM(3) 


No  <N<NT  -  V. 


N  =  N  +  1 
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From  255  or  from  after  271 


Return  to  1  for  next  case 
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C.5  PROGRAM  LISTING 


08300  C  ANALYTIC  TEMPERATURE  PROFILE  22  MARCH  1964 

003HO  20  FORMAT(2F10.5,2XI3,2X|2.2X14) 

08376  DIMENSION  T(24),  R{24),  A(7),  RJ(24),  RLUM{3),  TAVG(3),  F(17) 

08376  22  FORMAT (F5.1) 

08398  23  FORMAT (  F5 .0 , 1 5 , 2F1 0  «  5 ) 

08436  24  FORMAT (El 4.5) 

08458  151  FORMAT (8X22HARCRAD  PROFILE  NUMBER  ,14) 

08554  152  F0RMAT(8x5HN0RD«,  I2,3X3HRD=>,f6.3,3X3HRE=,f6.3) 

08678  153  FORMAT(8X5HRTOL«»,F8.6,3X5HTTOL~,F8.6) 

08776  275  FORMAT(9X40HN  TEMP  RADIUS  Rl  RO) 

08906  281  F0RMAT{7X13,3XF6.2,4XF7.4,3XF8.5(3XF8.5) 

09024  300  FORMAT(5XI2,3XF8.5,3XF9.5) 

09100  530  FORMAT(5X8hTMAX=B0=,F11-.5,3X3HB2=,F11.5) 

09194  531  FORMAT  (5X3mE0*«,  F1 1 ,5,3X3HE1b,F1 1 .5,3X3HE2=,F1  1.5) 

09308  532  FORMAT(5X3HTO=,F9.5,3X4HTDP-,F10.5,3X3UTE=,F9.5,3X4HTEP=*,F10.5) 

09458  1  READ  20,Rn,REp  NT,  WORD,  PRNUM 

09530  NCST=NORD+1 

09578  :i=1 

09616  10  READ  24,  A(N) 

09662  f|-N+l 

09710  I  F(IICST-IJ)  15,  10,  10 

09778  15  N-1 

09814  16  READ  22,  T(N) 

09862  N-N+l 

09910  IF(NT-N)  25,  16,  16 

09978  25  READ  23,  DELT,  IMAX,  RTOL,  TTOL 

TO038  5025  L=1 
T0074  RB=RD 

TP110  _ 2iRPOV/=1. 
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T0 146 
TO  1 82 
T021B 
T0254 
T0266 
T0314 
T041O 
T0458 
T057O 
T3646 
T0682 
T07IO 
T075^ 
T0790 
T0798 
T0834 
TH870 
Tig  54 
T  1  1!  1  4 
Tl  194 
T  7278 
T1338 
T 1 342 
T1410 
T1443 
Tl  466 
T1502 
T1550 
T 1 6 1 0 
T1646 
T1706 
T17S4 
T 1046 


TNR=A(1) 

TPMR=0. 

COEF=fl. 

DO  27  K-2,  HOST 
COEF  =C0£F+1. 

TPNR=TPNR+A(K)*RPOW*COEF 

RPOW*»RPOW*RB 

27  TI!R“TNR+RPOW*A  ( K) 

IF(L-l)  1000,  28,  29 

28  TD=TNR 
TDP=TPIIR 
RB-RE 
L=2 

GO  TO  26 

29  TE*=1NR 
TET  =TPNR 
[!0=TO-T3P*RD/  2. 

G2=TDP/2./RD 

£2=(  (TEP*(RE-1 .) )-(TE-T( 1 )) )/ (RE-1 ,  )**2 
E 1  “T  E  P-  2«*  E  2*R  E 
_  E0=T( 1 J-E2-E1 

30  PRINT  151,  PR HUM 

□I  PRINT  152,  IIORD,  RD,  RE 

32  PRINT  153,  RTOL,  TTOL 

33  IF(SENSE  SWITCH  3}  36,  35 

35  PRINT  530,  00,  B2 
PRINT  531,  E0,  El,  E2 
PRINT  532,  TO,  TOP,  TE,  TEP 

36  H*»  1 
R(N)-1- 

40  »=!.•+! 

44  IF(TE-T(r0)  60,  45,  45 

65  R( II)=-E1/ 2,/E2+SQRTF(  (El  /  2,/E2J**2+(T (N)-E0)/E2) 
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12098 

GO  TO  40 

T2106 

60  R( N)“RE 

T  2 1 66 

61  1 F(TD-T(N) )  100,  81,  62 

T2258 

62  TNR=A(1)-T(N) 

T2330 

TPNR=G. 

T2366 

RPOW-1. 

T2402 

COEF=0. 

T2438 

DO  63  K-2,  NCST 

T2450 

C0EF=COEF+1 . 

T2498 

TPNR=TPNR+A ( K)*RP0W*C0EF 

T2594 

RP0W“RP0W*R(H) 

T2666 

63  TNR=TNR+RPOW*A(K) 

T27B6 

RD!F=TMR/TPNR 

T2334 

R(N)“R(H)— RD1F 

T2930 

RDlF-SQRTF(RDIF*RDIF> 

T2990  i 

70  IF(RDIF-RTOL)  90,  90,  62 

T3058  | 

81  R(N)~RD 

T3118 

11=11+1 

T3166 

GO  TO  100 

T3174 

90  N-N+1 

T3222 

R(N)»R(H-1J 

T3308 

GO  TO  61 

T3314 

100  1 F(B0— T(N) )  101,  101,  110 

T3406 

101  MT-H 

T3442 

T ( NT) =B0 

T3502 

R(N)“0. 

T3562 

GO  TO  190 

13570 

110  R(N)“SQRTF((B0-T(N))/B2) 

T3690 

N“N+1 

T3738 

GO  TO  100 

T3746 

190 1)0  200  1=1,  17 

T3758 

200  F(  l  )*»0. 

T3854 

RJ(  1 ) =R{ 1 ) 

Statements  added  for  arbitrary  profiles. 


Eq.  3.  Z 

P  = NORD 

Q  =  SQRTF{A(NCST)*A(NCST)) 

DO  110,  N=l,  NT 
l100  R(N)  =  (l.  -  (T(N)-Z.  )/Q)*5^l,/P) 
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T3890 
73926 
73962 
74022 
74082 
74130 
T4214 
T4274 
T4286 
T4354 
T4422 
T4458 
T4494 
745  06 
74554 
74674 
74758 
74766 
74086 
74894 
75062 
75146 
T521O 
75338 
75374 
75410 
T5482 
75530 
75598 
T5754 
75934 
75954 
7605(5 


N-2 

202  L-0 

RULIM-RJ(H-l) 

RLLIM(  1)=»R(N) 

205  L-L+1 

DRAD=(RUL I M-RLL I M( L) ) / DELT 
RAD=RLLIM(L) 

DO  220  1-1,  (MAX 

214  I F(RAD-RD)  217,  217,  215 

215  I F(RAD-RE)  216,  218,  218 

216  FP0W=A( 1 ) 

RPOW=1. 

DO  2161  K=2,  NCST 
RPOW-RPOW*RAO 
2161  FPOVf=FPOW+RPOV/*A{K) 

F{  1 )  =FPOV(**4*RAD 
GO  TO  220 

217  F ( 1 ) = ( ( D0+8  2*RAD**2)**4)*RAD 
GO  TO  220 

218  F(l )=((E0+E1*RAD+E2*RAD**2)**4)*RAD 
220  RAD  =RAD+DRAD 

SUM1=F( 1 )+F( IMAX} 

SUf‘2=F(  2)+F(4)+F(6)+F{8)+F(10)+F(  12)+F(14)+F(16) 
SUM3-0. 

I  “3 

226  SUM3“SUM3+F(  I ) 

1  =  1+2 

I F( 1-1  MAX)  226,  227,  227 

227  FI UTG=.66666667*DRAD* (SUH1+4.*SUM2+2.*SUM3) 
TAVG(L)=(F! NTG/ {RULIM**2-RLLI M(L)**2)}**0.25 
I F(SENSE  SWITCH  2)  228,  229 

228  PRINT  300, L,  RLLIM(L),  TAVG(L) 

229  I F(L-2)  230,  240,  250 
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T6118  230  RLUM(2)=R(N+1) 

T61 78  GO  TO  205 

T61B6  240  CWTAVG(L-1}-TAVG(L))/(RLLIK(L-1)-RLUM(L)) 

76378  C0"TAVG(L-1 )-Ct*RLLIM(L-1) 

76498  RLL IM(3)=(T(N)-C0)/C1 

T6582  I  F(RLLIM(3))  241,  205,  205 

T6638  241  RLUM(3)»0. 

76674  GO  TO  205 

T6682  250  TDIF-TAVG(L)-T(M) 

T6778  253  IF(RLLIM{3))  1000,  254,  258 

T6B34  254  JF(TDIF)  255,  258,  258 

76890  255  RJ(N)=0. 

T6950  IF(HT-IM)  280,  280,  256 

17030  256  N-N+1 

T7078  GO  TO  255 

T7086  258  TDI FA=SQRTF(TDI F*TOI F) 

T7146  I F(TDI FA-TTOL)  270,  270,  260 

T7214  260  00265  K»1,2 

T7226  TAVG(K)-TAVG(K+1} 

T7310  265  RLUM(K)-RLLIM(K+1) 

T7430  L-L-l 

T7478  GO  TO  240 

T7486  270  IF(RLLIM(3)-R(N+1))  272,  271,  271 

T7578  271  RJ{N)=.RLLIH(3) 

T7638  I F( NT— H-1 }  280,  280,  273 

T77 1 8  272  RJ(N}»  0. 

T7778  N-N+1 

T7826  GO  TO  270 

T7834  273  M-M+1 

77882  GO  TO  202 

77890  280  RJ(M>0. 

T7950  PRINT  275 

77974  N-1 
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18(110  PRINT  281,  14,  T  ( N ) ,  R(N) 

T8106  N=2 

181^2  282  PRINT  281,  tJ,  T(N) ,  R(N) ,  RJ(N),  RJ(N-I) 

T8310  N-N+1 

T83S8  I  F{NT“f4)  1  ,  282,  282 

T8426  1000  STOP 

T8U74  END 

=ROG  SW  1  ON  FOR  SYMBOL  TABLE,  PUSH  START 

51/  1  OFF  TO  IGNORE  SUBROUTINES,  PUSH  START 

PROCESSING  COMPLETE 
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APPENDIX  D 

CURVE  FIT  FOR  TEMPERATURE  PROFILE  DATA 
D.l  ARCRAD  CURVE  FIT  PROGRAM  DESCRIPTION 

Reduction  of  temperature  profile  data  to  analytic  form  was 
accomplished  by  making  use  of  a  modification  of  a  standard  computer 
program  by'W.  R.  Greaves  (Ref.  19)  entitled  POLYNOMIAL  CURVE 
FITTING,  NO  CONSTRAINTS. 

The  following  changes  were  made  in  Greaves'  program. 

1.  Format  statements  were  added  to  change  from  "No 
Format"  FORTRAN  to  FORTRAN  I. 

2.  The  maximum  order  polynomial  was  reduced  from 
50  to  7. 

3.  Provision  for  weighting  data  was  eliminated. 

4.  Control  statements  were  rewritten  to  the  form  shown 
in  the  accompanying  block  diagram. 

The  ARCRAD  CURVE  FIT  program  fits  a  straight  line  and 
successively  higher  order  curves  to  the  input  data  on  a  "least  squares" 
basis.  For  each  order  curve  the  polynomial  constants  and  the  standard 
error  are  printed.  As  the  order  of  the  fitted  curve  increases  the 
standard  error  at  first  decreases,  then  levels  off  and  either  increases 
or  oscillates.  A  curve  which  is  the  best  compromise  between  minimum 
order  and  minimum  standard  error  is  selected  as  most  representative 
of  the  data.  The  program  automatically  causes  a  card  to  be  punched 
for  each  constant  in  the  polynomial  for  use  as  input  to  the  ANALYTIC 
TEMPERATURE  PROFILE  program.  Finally,  for  the  order  selected, 
the  value  of  the  polynomial  at  each  input  data  point  is  printed,  along 
with  the  input  data  value  and  the  difference. 

The  number  of  data  points  must  be  at  least  2  greater  than  the 
order  of  the  maximum  order  polynomial.  The  program  can  handle  up 
to  50  sets  of  data  and  computes  polynomials  up  to  seventh  order  but 
these  limits  can  easily  be  changed  by  appropriate  changes  in  the 
DIMENSION  statement. 
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D.2  SENSE  SWITCH  SETTINGS 

Flexibility  in  program  utilization  is  made  possible  by  the  Sense 
Switches.  Switch  1  OFF  causes  automatic  selection  of  the  3rd  order 
polynomial.  Switch  1  ON  causes  the  computer  to  pause  after  printing 
the  constants  and  the  standard  error  for  the  preselected  maximum 
order.  The  desired  order  polynomial  can  then  be  selected  and 
inserted  in  the  machine  by  on-line  typewriter  (ACCEPT  statement). 

Switch  2  ON  eliminates  automatic  card  punching  of  the  poly¬ 
nomial  constants. 

D.3  ARCRAD  CURVE  FIT  DEFINITION  SYMBOLS 


AM  (I,  J) 

Matrix  of  coefficients  of  normal  equations 

DIF 

Difference  between  Y(I)  and  YFLEX(I) 

LAST 

Highest  order  polynomial  to  be  fitted 

NDATA 

Number  of  data  points 

PRNUM 

Serial  number  of  temperature  profile 

SUMX(J) 

N  J-l 

£  (Xif  1 

i=l 

SUMY(J) 

N  J-l 

S  (Yi)  (Xi)J 

i=l 

SI 

Same  as  YFLEX(I) 

I  N  2 

.1  (Yi  -  YFLEXi) 
Standard  error, V  NDATA-NORD-l 

S2 

S3 

Same  as  DIF 

X(I) 

th 

i  value  of  independent  variable  (radius) 

Y(I) 

i**1  value  of  dependent  variable  (temperature) 

YFLEX(I) 

Value  of  the  polynomial  at  X(I) 
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D.4  BLOCK  DIAGRAM 


2001 
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from  2093 


(hich  time  tKrnutpn  secor 
he  program? 


Print:  NORD, 
Sigma,  Curve 
constants 


T 


2161 

No  <r*NORP  =  LasT 


Print  NORD,  Sigma 
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D.5  PROGRAM  LISTING 


08300  C  ARCRAD  CURVE  FIT  10  MARCH  1964 

08300  DIMENSION  X(52),  Y(52),  AM(9,9),  SUMX{15),  SUMY(17) 

08300  2501  F0RMAT(I2, 2X13, 2X14) 

08354  2502  FORMAT ( F6 , 2 ,  2X  E 1 4.7} 

38392  2503  FORMAT() 

38410  2504  FORMAT (5X5HN0RD-,  I2,3X6HS1GMA-,F7.4) 

38504  2505  FORMAT (5X1 2, 3XF1 1.5) 

38562  2506  FORMAT (8X11HCASE  NUMBER, l 4,3x6HNDATA=, I  3) 

08674  2507  FORMAT(5XF7.4,3XF6.2,3XF8.4,3XF7.4) 

08768  2508  FORMAT (5X47HSYJ1  OFF  LAST-3,  SW1  ON  ACCEPT,  LAST,  PUSH  START) 

08904  2509  F0RMAT(7X34HX(I)  Y(l)  YFLEX(I)  DIF) 

09018  2510  F0RMATCE14.5) 

09040  2001  READ  2501,  LAST,  NDATA,  PRNUM 

09088  PRINT  2506,  PRNIJM,  NDATA 

09124  NRPT—0 

09148  1-0 

09172  2020  1=1+1 

09208  READ  2502,  X(1),  Y(l) 

09292  I F( I -NDATA)  2020,  2070,  2070 

39360  2070  SUMX( 1 )-fl, 

39384  SUMX(2)=0. 

39408  SUMX(3)=0. 

09432  SUMY( 1 )=0. 
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09456 

09480 

09492 

09528 

09588 

09684 

09744 

09876 

09900 

09936 

09972 

09904 

09996 

T0044 

T01J38 

T0332 

T0344 

13440 

10476 

T0408 

T060O 

T0692 

70704 

70716 

T1088 

T 1  TOO 

Til  12 

T1364- 

T1388 

T1400 

71424 

T 1 508 

Tl52fl 


$UMY{2)=0. 

DO  2093  1=1,  HDATA 
suMX(f)=suMx(i)+r; 
SUMX(2)=SUMX(2)+X(I) 
SUMX(3)=SUMX(3)+X( I )*X(f ) 

SUMY( 1)=SUMY(1)+Y(I ) 

2093  SUMY(2)»SUMY(2)+X(1 )*Y(I ) 

N0RD=1 
2093  L-NORD+1 
KK-L+1 

DO  2101  |=1 ,  L 
DO  2100  J=1,  L 
I K-J-1+1 

2100  AM(  I , J)-SUMX(  1 K) 

2101  AI4(  I  ,KK)=SU11Y(I  ) 

DO  2140  1=1,  L 
AM(KK,I)  — 1. 

KKK-I+I 

DO  2110  J=KKK,  KK 
2110  AM(KK,J)=0, 

C-1./AM(1,I) 

DO  2120  11=2,  KK 
DO  2120  J=KKK,KK 

2120  AM( M , J)=AH( I  I , J)-AM( 1 , J)*AM( I  I , I }*C 
DO  2140  11=1,  L 
DO  2140  J=KKK,KK 
2140  AM(I I , J)=AM( I  1  +  1, J) 

S2-0. 

DO  2160  J=1,  IJDATA 
S1-0. 

S 1=  S1+  AM(1,  KK) 

DO  2150  1=1,  HORD 
2150  S 1-S 1+AM( 1+1 ,KK)*X( J)**l 
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7 1700 
T 1880 
T1928 
71976 

T2032 

72100 

72124 

T2160 

72184 

72220 

72256 

72352 

72420 

72488 

72508 

72544 

72568 

72588 

72612 

72648 

72732 

72800 

72824 

72860 

72884 

72956 

72968 

T3148 

73208 

73316 

73384 

73420 

T3456 


2160  S2=S2+(Sl-Y( J) )*($1-Y(J)} 
B-NDATA-L 

S2=(S2/B)**';5 
IF(NRPT)  2163,  2163,  2161 

2161  IF(NORD-UST)  2171,  2167,  2167 

2163  PRINT  2503 

PRINT  2504,  WORD,  S2 
1-0 

2164  1-1+1 
J- 1-1 

PRINT  2505,  J,  AM( I ,KK) 

IF(I-L)  2164,  2165,  2165 

2165  IF  (NORD-LAST)  217.1,  2166,  2166 

2166  I F(SENSE  SWITCH  1}  2173,  2167 

2167  PRINT  2504,  NORD,  S2 
PRINT  2509 

I F( SENSE  SWITCH  2)  2203,  2201 
220!  1=0 
2202  1-1+1 

PUNCH  2510,  AM( I ,  KK) 

IF(I-L)  2202,  2203,  2203 

2203  1=0 

2204  1=1+1 

S  1=0.' 

$1=  AM(1,KK) 

DO  2168  J-1,  NORD 

2168  S 1°S 1+AM( J+1 ,KK)*X{ I )**J 

S3=Y( I }-Sl 

PRINT  2507,  X(l),  Y( I ) ,  SI,  S3 
IF(I-NDATA)  2204,  2001,  2001 
2171  NORD-NORO+1 
J=2*N0RD 
SUMX( J)=0. 
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T3504  SUMX(J+1)-0. 

T3552  SUMY(N0RD+1)»fl. 

T3600  00  2172  1-1,  HDATA 

T3612  SUMX(J)-SUMX{J)+X(I)**(J-1) 

T3768  SUMX(J+1  }s»SUMX(J+1)+X(l  )**J 

T38B8  2t72  Sl/MY(HORW-t)»SUMY(NORD+1)+YO)*X(l)**NORD 

T40B0  GO  TO  2093 

T4088  2173  NRPT-NRPT+1 

T4124  PRINT  2508 

T4148  PAUSE 

T4160  IF  (SENSE  SWITCH  1)  2175,  2176 

T4180  2175  ACCEPT  2501,  LAST 

T4204  GO  TO  2070 

T4212  2176  LAST=3 

T4236  GO  TO  2070 

14244  END 

PROG  SW  1  ONFOR  SYMBOL  TABLE,  PUSH  START 
SW  1  OFF  TO  IGNORE  SUBROUTINES,  PUSH  START 

PROCESSING  COMPLETE 
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APPENDIX  E 

MEASUREMENT  OF  TEMPERATURE  PROFILES 


Conditions  required  for  application  of  the  intensity  ratio  tech¬ 
nique  are:  (1)  steady  symmetric  Bource,  (2)  Boltzmann  distribution 
of  the  atomic  energy  levels  (i.  e. ,  local  thermal  equilibrium  of  the 
plasma  particles)  and  (3)  optically  thin  plasma  at  the  frequency  of  the 
emission  line.  The  first  condition  was  checked  photographically  with 
a  Fastax  camera  operating  at  3000  to  5000  frames  per  second.  In  this 
speed  range,  which  was  of  the  same  order  as  the  spectrograph 
exposures  of  1/1000  sec,  the  Bource  appeared  steady  and  uniform. 

The  existence  of  local  thermal  equilibrium  is  discussed  in  Section  2.  1. 
Optical  thinness  is  one  of  the  factors  governing  the  choice  of  suitable 
lines.  At  one  atmosphere  severalglines  can  be  used,  but  at  2  atmos¬ 
pheres  only  the  oxygen  pair  7774  A  -  7947  A  is  suitable  and  it  iB  near 
marginal.  Relative  intensity  at  a  given  wavelength  is  related  to 
transition  probability  A,  atom  number  density  N  and  photon  energy 
hP  by 


r\l  A1N1(hl/)1 

XX2  =  A2N2(hy)2 


(E-l) 


For  spectral  lines  of  the  same  species, 
tribution  of  energy  states 


TX1  =  gl  Aixi 
TXZ  g2  A2  ^2 


assuming  a  Boltzmann  dis- 


(E-2) 


where  g  is  the  statistical  weight  of  an  energy  state  and  E  is  the 
excitation  energy  of  the  upper  state  energy  level.  Solving  for  T 


(Ej  -  Eg)  (1/k) 

ln*glAl 


(E-3) 
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The  excitation  energies  and  transition  probabilities  are  known,  so 
spectrographic  measurement  of  the  intensity  ratios  establishes  an 
"average"  temperature  of  the  column.  To  determine  the  temperature 
distribution  the  true  intensity  distribution  must  be  calculated  from 
the  "average"  observed  intensity  using  Abel's  integral  equation 


I(x) 


f 

J  T  —X 


n(r)  rdx 

yj  r^  - 


(E-4) 


where  I(x)  is  observed  intensity  over  radius  x  and  n(r)  is  the  number 
density  of  emitters  in  particular  radial  bands  of  radius  r.  This 
function  has  been  tabulated  by  Pearce  (Ref,  20). 

Temperature  level  can  be  established  because  the  spectral  line 
intensity  goes  through  a  maximum  at  a  unique  temperature  for  a 
given  pressure.  This  maximum,  the  so  called  peaking  function, 
results  from  the  counter -balancing  effects  of  increased  excitation  of 
the  upper  states,  which  increases  the  number  of  emitters  per  unit 
volume,  and  increased  ionization  which  reduces  the  number  of 
emitters  per  unit  volume.  The  technique  has  been  described  in  detail0 
by  Sadjiar^tRef,  21  ),  For  this  study  the  oxygen  atom  pair,  X  =  7947  A 
and  7774  A  were  used well  as  the  nitrogen  atom /ion  pair, 

\  =  4137  A  and  3995  A. 
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